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ABSTRACT 


Linear least squares estimation and regression analyses continue to 
play a major role in orbit determination and related areas. In this report 
we document a library of FORTRAN subroutines that have been developed to 
facilitate analyses of a variety of estimation problems. Our purpose is to 
present an easy to use, multi-purpose set of algorithms that are reasonably 
efficient and which use a minimal amount of computer storage. Subroutine 
inputs, outputs, usage and listings are given, along with examples of how 
these routines can be used. The following outline indicates the scope of 
this report: Section I, introduction with reference to background material 

Section II, examples and applications; Section III, a subroutine directory 
summary; Section IV, the subroutine directory user description with input, 
output and usage explained; and Section V, subroutine FORTRAN listings. 

The routines are compact and efficient and are far superior to the normal 
equation and Kalman filter data processing algorithms that are often used 
for least squares analyses. 
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I. 


Introduction 


Techniques related to least squares parameter estimation play a 
prominent role in orbit determination and related analyses. Numerical 
and algorithmic aspects of least squares computation are documented 
in the excellent reference work by Lawson and Hanson, Ref. [1]. Their 
algorithms, available from the JPL subroutine library, Ref. [2], are 
very reliable and general. Experience has, however, shown that in 
reasonably well posed problems one can streamline the least squares 
algorithm codes and reduce both storage and computer times. In this 
report, we document a collection of subroutines most of which we have 
written that can be used to solve a variety of parameter estimation 
problems . 

The algorithms for the most part involve triangular and/or 
symmetric matrices and to reduce storage requirements these are stored 
in vector form, e.g. , an upper triangular matrix U is written as 


U 11 U 12 U 13 

U 14 


0(1) 

U(2) 

U(4) 

U<7) 

U 22 U 23 

U 24 , 

etc. 

. 


U(3) 

U(5) 

etc. 

A 

D 34 ' 


A 


U (6) 

U(9) 

(J 

U 44 


(J 



U(10) 

l V/ 






- 


Thus, the element from row i and column j of U, i < j, is stored in 
vector component j(j-l)/2 + i. We hasten to point out that the engineer, 
with few exceptions, need have no direct contact with the vector sub- 
scripting. By this we mean that the vector subscript related operations 
are internal to the subroutines, vector arrays transmitted from one 
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subroutine to another are compatible,, and vector arrays displayed 
using the print subroutine TWOMAT appear in a triangular matrix format. 

Aside : The most notable exception is that matrix problems are generally 

formulated using doubly subscripted arrays. Transforming a double 

subscripted symmetric or upper triangular matrix A(-j-) to a vector 

stored form, U(-) is quite simply accomplished in FORTRAN via 
IJ = 0 

DO 1 J = 1,N 
DO 1 I = 1,J 
IJ = IJ+1 
1 U(IJ) = A(I, J) 

Similarly, transforming an initial vector D(*) of diagonal positions of 
a vector stored form, U('), is accomplished using 



JJ = 0 


JJ = N*(N+l)/2 


DO 1 J = 1,N 

or 

DO 1 J = N,l,-1_ 


JJ = JJ+J 


✓"-s 

v/ 

1 ! 

/-— \ 
*”3 

s— ' 

1 

U(JJ) = D (J) 


i jj = jj-j 


The conversion on the right has the modest advantage that D and U can 
share common storage (i.e., U can overwrite D) . These conversions 
are too brief to be efficiently used as subroutines. It seems that when 
such conversions are needed one can readily include them as in-line code. 

End of Aside 

This package of subroutines is designed, in the main, for the analysis 
of parameter estimation problems. One can, however, use it to solve problems 
that involve process noise and to map (time propagate) covariance or infor- 
mation matrix factors. In the case of mapping the storage savings associated 
with the use of vector stored triangular matrices is, to some extent, lost. 
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Mathematical background regarding Householder orthogonal trans- 
formations for least squares analyses and U-D matrix factorization 
for covariance matrix analyses are discussed in references [1] and [3]. 

Our plan is to illustrate, in Section II, with examples, how one can 
use the basic algorithms and matrix manipulation to solve a variety 
of important problems. The subroutines which comprise our estimation 
subroutine package are described in Section III, and detailed input/ 
output descriptions are presented in Section IV. 

Section V contains FORTRAN listings of the subroutines. There are 
several reasons for including such listings. Making these listings 
available to the engineer analyst allows him to assess algorithm 
complexity for himself; and to appreciate the simplicity of the 
routines he tends otherwise to use as a black box. The routines use 
only FORTRAN IV and are therefore reasonably portable (except possibly 
for routines which involve alphanumeric inputs) . When estimation 'problems 
arise to which our package does not directly apply (or which can be made 
to apply by an awkward concatenation of the routines) one may be able to 
modify the codes and widen still further the class of problems that can be 
efficiently solved. 
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II. APPLICATIONS AND EXAMPLES 


Our purpose in this section is to illustrate, with a number of examples, 
some of the problems that can be solved using this ESP. The examples, in 
addition, serve to catalogue certain estimation techniques that are quite 
useful. 


To begin, let us catalogue the subroutines that comprise the ESP: 


1) A2A1 

2) COMBO 

3) COVRHO 

4) C0V2RI 

5) C0V2UD 

6) C2C 

7) INF2R 

8) HHPOST 

9) PERMUT 

10) PHIU 

11) RA 

12) RANK1 

13) RCOLRD 

14) RINCON 

15) RI2C0V 

16) R2A 

17) R2RA 

18) RUDR 

19) SFU 

20) TDHHT 

21) THH 

22) TTHH 

23) TWOMAT 


(A to A one) 
-(combo) 

(cov rho) 
(cov to RI) 
(cov to U-D) 
(C to C) 

(inf to R) 
(HH POST) 
(permut) 
(PHI*U) 

(R*A) 

(rank 1) 

(R colored) 
(r in- con) 

(Rl to cov) 
(R to A) 

(R to RA) 

(rudder) 

(S F U) 

(TDHHT) 
(T H H) 

(TTHH) 

(two mat) 


Matrix A to matrix A1 

Combine R and A namelists 

Covariance to correlation matrix, RHO 

Covariance to R inverse 

Covariance to U-D covariance factors 

Permute the rows and columns of matrix C 

Information matrix to (triangular) R factor 

Householder triangularization by post multiplication 

Permute the columns of matrix A 

Multiplies a rectangular PHI matrix by the vector 
stored U matrix that has implicitly defined unit 
diagonal entries . 

R(upper triangular, vector stored) *A (rectangular) 

Updated U-D factors of a rank-1 modified matrix 

(SRIF)R colored noise time-update 

R inverse along with a condition number bounding 
estimate 

R inverse to covariance 

Triangular R to (rectangular stored) matrix A 

Transfer to triangular block of (vector stored) R 
to a triangular (vector stored) RA 

(SRIF)R to U-D covariance factors, or vice-versa 

Sparse F matrix * vector stored U matrix with 
implicitly defined unit diagonal entries 

Two dimensional Householder matrix triangularization 

Triangular vector stored Householder data processing 
algorithm 

Orthogonal triangularization of two triangular 
matrices 

Two dimensional labeled display of a vector stored 
triangular matrix 
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24) 

TZERO 

(T zero) 

Zero a horizontal segment of a vector 
stored triangular matrix 

25) 

UDCOL 

(U-D colored) 

U-D covariance factor colored noise update 

26) 

UDMEAS 

(U-D measurement) 

U-D covariance factor measurement update 

27) 

UD2C0V 

(U-D to cov) 

U-D factors to covariance 

28) 

UD2SIG 

(U-D to sig) 

U~D factors to sigmas 

29) 

UTINV 

(U T inverse) 

Upper triangular matrix inverse 

30) 

UTIROW 


Upper triangular inverse, inverting only 
the upper rows 

31) 

WGS 

(W G-S) 

U-D covariance factorization using a weighted 
Gram-Schmidt reduction 


These routines are described in succeedingly more detail in sections III, 
IV, and V. The examples to follow are chosen to demonstrate how these 
various subroutines can be used to solve orbit determination and other 
parameter estimation problems. It is important to keep in mind that these 
examples are not by any means all inclusive, and that this package of 
subroutines has a wide scope of applicability . 

II. 1 Simple Least Squares 

Given data in the form of an overdetermined system of linear 
equations one may want a) the least squares solution; b) the estimate 
error covariance, assuming that the data has normalized errors; and 
c) the sum of squares of the residuals. The solution to this problem, 
using the ESP can be symbolically depicted as 

• [A:a] JSL[R:z], e 

Remarks : The array [A:z] corresponds to the equations Ax = z-v, veN(0,I); 

A A A ^ A 

the array [R:z] corresponds to the triangular data equation Rx = z-v, 
vsN(0,I) and e = ||z-Ax|| 




[R:z] 


utinv r ;-i % 

►[R : x] 


Remark : 



A 

Z 
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One may be concerned with the integrity of the computed inyerse 
and the estimate. If one uses subroutine RINCON instead of UTINy then 

m addition one obtains an estimate (lower and upper bounds) for the 

condition number R. If this condition number estimate is large the 

computed inverse and estimate are to he regarded with suspicion. By 

large, we mean considerable with respect to the machine accuracy (viz. 

on an 18 decimal digit machine numbers larger than 10^) . Note that the 

condition number estimate is obtained with negligible additional compu- 


tation and storage. 


'"-1 
ER 1 


RI2C0V [C] 


Remarks : C = R R = estimate error covariance. Some computation can 

be avoided in RI2C0V if only some (or all) of the standard deviations 
are wanted. 

II. 2 Least Squares With A Priori 

If a priori information is given, it can be included as additional 
equations (in the A array) or used to initialize the R array in subroutine 
THH (see the subroutine argument description given in section IV) . One is 
sometimes interested in seeing how the estimate and/or the formal 
statistics change corresponding to the use of different a priori 
conditions. In this case one should compute [R;z] as in case II. 1, and 
then include the a priori [R^iz^] using either subroutine THH, or 
subroutine TTHH when the a priori is diagonal or triangular, e.g., 

A A 

[R:z] 


[R:*n 


[R 


o o 


TTHH 


The new result overwrites the old. 
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It is often good practice to process the data and form [R:z] before 
including the a priori effects. When this is done one can analyze 

the effect of different a priori, [R q : z q ] without reprocessing the data. 

If a priori is given in the form of an information matrix. A, 

(as for example would be the case if the problem is being initialized 

k 

with data processed using normal equation data accumulation ) then one 
can obtain R q from A using INF2R; 


A 


INF2R 


o 

f T 1 ^TTl 

If there were a normal equation estimate term, z = A b, then z =R z. 

o o 

II. 3 Batch Sequential Data Processing 

Prime reasons for batch sequential data processing are that many 
problems are too large to fit in core, are too expensive in terms of core 
cost, and for certain problems it is desirable to be able to incorporate 
new data as it becomes available. Subroutines THH and UDMEAS are specially 
designed for this kind of problem. Both of these subroutines overwrite 
the a priori with the result which then acts as a priori for the next 
batch of data. If the data is stored on a file or tape as A^ z^, A 2 , z 2 ,... 
then the sequential process can he represented as follows : 

SRIF Processing* * 

a) Initialize [R:z] with a priori values or zero 

b) Read the next [A:z] from the file 


* T T T 

i.e., solving Ax = b-v with normal equations, A Ax q = A b; A = A A 

is the information matrix. 

The acronym SRIF represents Square Root Information Filter. The SRIF is 
discussed at length in the book by Bierman, ref. .[3]. 
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d) If there is more data go back .to b) 

e) Compute estimates and/or covariances using UTINV and RI2C0V 
(as in example II. 1) 


’EM)** Processing 

A A A 

a') Initialize [U-D:x] with a priori U-D_ covariance factors and the 
initial estimate 

b^) Read the next [A:z] scalar measurement from the file 


o 

d') 

O 


[U-D:x] 


UDMEAS 


[U-D: x] 


[A: z] } 

If there is more data go back to b') 

Compute standard deviations or covariances using UD2SIG or 


UD2C0V. 


Note that subroutine THH is best (most efficiently) used with 
data batches of substantial size (say 5 or more) and that UDMEAS processes 
measurement vectors one component at a time. If the dimension of the 
state is small the cost of using either method is generally negligible. 

The UDMEAS subroutine is best used in problems where estimates are 
wanted with great frequency or where one wishes to monitor the effects 
of each update. In a given application one might choose to process 
data in batches for a while and during critical periods it may be 


The new result overwrites the old . 

** 

U-D processing is a numerically stable algorithmic formulation of the Kalman 
filter measurement update algorithm, cf reference [3]. The estimate error 
covariance is used in its UDU T factored form, where U is unit upper triangular 
and D is diagonal. 
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desirable to monitor the updating process on a point by point basis. 

In cases such as this, one may use RUDR to convert a SRIF array to U-D 
form or vice-versa. 

Remarks : Another case where an R to U-D conversion can be useful occurs 

in large order problems (with say 100 or more parameters) where after 
data has been SRIF processed one wants to examine estimate and/or 
covariance sensitivity to the a priori variances of only a few of the 
variables. Here it may be more convenient to update using the UDMEAS 
subroutine. 

II. 4 Reduced State Estimates and/or Covariances From a SRIF Array 

Suppose, for example, that data has been processed and that we have a 
A A 

triangular SRIF array [R:z] corresponding to the 14 parameter names, a^, a^., 
a , x, y, z, v , v , v , GM, CU41, L041, CU43, L043 (constant spacecraft 

y» * J 9 9 x y z 

accelerations, position- and velocity, target body gravitational constant, 
and spin axis and longitude station location errors). 

Let us ask first what would the computed error covariance be of 
a model containing only the first 10 variables, i.e. , by ignoring the 
effect of the station location errors. One would apply UTINV and RI2C0V 
just as in example II. 1, except here we would use N (the dimension of 
the filter ) = 10, instead of N=14. 

Next, suppose that we want the solution and associated covariance 
of the model without the 3 acceleration errors. One ESP solution is to 
use 
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errors a , a , or a . 

r x’ y 

Note how triangularizing the rearranged R matrix produces the 
desired lower dimensional SRIF array; and this is the same result one 
would obtain if the original data had been fit using the 11 state model. 

As the last subcase of this example suppose that one is only 
interested in the SRIF array corresponding to the position and velocity 
variables. The difference between this example and the one above is 
that here we want to include the effects due to the other variables. 


z is often given the label RHS (right hand side) 
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One might want this sub-array to combine with a position-velocity SRIF 


array obtained from, say, optical data. 


A A 

• [Rjz] 


R2RA 


INPUT NAMES: 

a _» a > a , x, y, z, v , v , v , GM 
y x y z 

CU41, L041, CU43, L043, RHS 


One method to use would be. 


[ r a :z a3 


OUTPUT NAMES: 


x, y, z, v x> v , v z , GM 
CU41, L041, CU43, L043, RHS 


Remark: The lower triangle starting with x is copied into R . 

il 


R2A 

[ R A :Z A^ :z A-* ( Reorderin g) 


NAMES: GM, CU41, L041, CU43, L043, 

x, y, z, v^, v , v z , RHS 

THH A A 

• [A jz^j *^ R a :Z A^ ( Trian 8 ularizij:1 8) 

* £ R a : z i} *"[R :z ] (Shifting array) 

xl A XX 

NAMES: x, y, z, v , v , v , RHS 

x y z 

Remark » The lower right triangle starting with x is copied into R^. 

We note that one could have elected to use COMBO in place of the first 
R2RA usage and R2A; this would have involved slightly more storage, but 
a lesser number of inputs. The sequence of operations is in this case. 


• [R:z] 


COMBO r . , 

*-[A:z] 


ORIGINAL NAMES 


DESIRED NAMES: 


x, y. 


z, v , v , v , RHS 
x y z 


Remark ; By using COMBO the columns of [R:z] , are ordered corresponding to 

t 

the names a , a , a , GM, CU41, L041, CU43, and L043, followed by -the 
r x y 

desired names list. 
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TTTH A A 

• [A: z [R:z] 

* A A 

Remark : The [R:z] array that is output from this procedure is 

equivalent but different from the [R:z] array that we began with. 


• [R:z] ^^ .[R :z ] 

xx 

Remark: As before, the lower right triangle starting with x is copied 


into R . 
x 


To delete the last k parameters from a SRIF array, it is not 
necessary to use subroutines R2A and THH. The first N - k = N columns 
of the array already correspond to a square root information matrix of 
the reduced system. If estimates are involved one can simply move the 
z column left using: 

R (N*(N + l)/2 + i) = R(N*(N +l)/2 + i) , i = l,...,k. 


Remark : We mention in passing that if one is only interested in estimates 

and/or covariances corresponding to the last k parameters then one can use 

R2RA to transform the lower right triangle of the SRIF array to an upper 

left triangle after which UTINV and RI2C0V can be applied. 

II. 5 Sensitivity, Perturbation, Computed Covariance and Consider 
Covariance Matrix Computation 

Suppose that one is given a SRIF array 



N 

R 

xy 



0 



\\ 


(II. 5a) 
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in which the N variables are to be considered. (One can, of course, using 

y 

subroutines R2A and THH reorder and retnangularize an arbitrarily arranged 
SRIF array so that a given set of variables fall at the end.) For various 
reasons one may choose to ignore the y variables- in the equation 


R x + R y = z 
x xy x 


v , v eN(0,I) 

A. A 


and take as the estimate x = R z . It then follows that 

c xx 


x - x 

c 




y - 
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(II. 5b) 


(II. 5 c) 


and from this one obtains 


Sen = 


3 (x-x c ) 
3y 



(II. 5d) 


(sensitivity of the estimate error to the unmodeled y parameters) 

Pert = Sen * Diag (o^(l) , . . . ,o^(K^)) (II. 5e) 

where a (1) , . . . ,a^(N^) are a priori y parameter uncertainties. 

(The perturbations are a measure of how much the estimate error could be 
expected to change due to the unmodeled y parameters.) 


P = R 1 R T + Sen P Sen T (II. 5f) 

con x x y 

T t 

*= + (Pert) (Pert) if P^ is diagonal 


where P is the estimate error covariance of the reduced model, 
c 

An easy way to compute P^, Pert and P is as follows: Use subroutine 
R2RA to place the y variable a priori [P^(0) : y^] into the lower right 


1 

Pert = Sen 



n 


The a priori estimate y 


of consider parameters is generally zero. 
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corner of (II. 5a), replacing R and z , i.e. : 


[R : z ] 


R2RA 


[P y 2 (0) 


^ 1 1 
y L 


R 

x 

0 


R 

xy 


P % (0) 

y 


z 

X 


A 

y„ 


Now apply subroutine UTIROW to this system (with a -1 set in the lower right 
comer* ) 


r 




— 


— 

R 

X 

! 

' P3 

1 

1 

z 

X 


R- 1 

X 

** 

Pert 

X 

c 

0 

P^(0) 

y 

A 

UTIROW r 

0 

p 3s (o) 

y 

A 

y o 


y o 

0 

0 

-1 


0 

0 

-1 

1— 


— 


— 




Note that the lower portion of the matrix is left unaltered, i.e., the purpose 

of UTIROW is to invert a triangular matrix, given that the -lower rows have 

already been inverted. From this array one can, using subroutine RI2C0V, 

get both P and P 

c con 

r — 1-, RI2C0V r-p -t „ , 

[R 1 LP j computed covariance 

x c 

RT2COV 

l\ ‘ Pert] [P con 3 consider covariance 

Suppose now that one is dealing with a U-D factored Kalman filter for- 
mulation. In this case estimate error sensitivities can be sequentially 

To have estimates from the triangular inversion routines one sets a -1 in the 
last column (below the right hand side) . 

** 

Strictly speaking this is not what we call the perturbation unless P v (0) is 
diagonal . ^ 
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Calculated as each scalar measurement (z = a^x + a^y + v) is processed. 

Sen. = Sen. , - K. (a T Sen. , + a T ) 

3 3-1 3 x 3-1 y' 

where Sen.. ^ is the sensitivity prior to processing this (j-th) measurement, 

■f 

and K. is the- Kalman gain vector. 

In this formulation one computes P ^ in a manner analogous to that des- 
cribed in- section II. 7; 

Let U x = U_. , = D.. (filter U-D factors) 

[s ] = Sen. (estimate error sensitivities) 

n y J 

then recursively compute 


U. -D 
k 


k 



KAHKl 


U k+l“ D k+l 


k = 1 , 


n 

y 


For the final U-D we have 


= U , D C °” = D 

3+1 n +1 ’ i+l n +1 

y y 


If P (0) = U D > instead of P (0) = Diag (a,,..., ) , then in the 

y yyy y l n y y * 

U— D recursion one should replace the Sen. columns by those of Sen.*U and 

J 3 y 

2 

0 \ should be replaced by the corresponding diagonal elements of D^. 

II . 6 Combining Various Data Sets 

In this example we collect several related problems involving data sets 
with different parameter lists. 

Suppose that the parameter namelist of the current data does not 
correspond to that of the a priori SRIF array. If the new data involves 
a permutation or a subset of the SRIF namelist, then an application of 

= g/a where g and a are quantities computed in subroutine UDMEAS. 
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subroutine PERMUT will create the desired data rearrangement. If the data 
involves parameters not present in the SRIF namelist then one could use 
subroutine R2A to modify the SRIF array to include the new names and then 
if necessary use PERMUT on the data, to rearrange it compatibly. 

Suppose now that two data sets are to be combined and that each 
contains par ame ters peculiar to it (and of course there are common para- 
meters) . For example let data set 1 contain names ABC and data set 2 
contain names DEB. One could handle such a problem by noting that the list 
ABCDE contains both name lists. Thus one could use subroutine PERMUT 
on each data set comparing it to the master list, ABCDE, and then the 
results could be combined using subroutine THH. An alternative automated 
method for handling this problem is to use subroutine COMBO with data 
set 1 (assuming it is in triangular form) and namelist 2 . The result 
would be data set 1 in double subscripted form and arranged to the name- 
list ACDEB (names A and C are peculiar to data set 1 and are put first) . 
Having determined the namelist one could apply subroutine PERMUT to data 
set 2 and give it a compatible namelist ordering. 

The process of increasing the namelist size to accommodate new 
variables can lead to problems with excessively long namelists, i.e., 
with high dimension. If it is known that a certain set of variables 
will not occur in future data sets then these variables can be eliminated 
and the problem dimension reduced. To eliminate a vector y from a SRIF 
array, first use subroutine R2A to put the y names first in the namelist; 
then use subroutine THH to retriangularize and finally use subroutine R2RA 
to put the y independent subarray in position for further use; viz. 
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m [ a ] 


R R z 
y yx y 


0 R z 
x x 


R2RA 


[R :z ] 
xx 


The rows [R^tR^^iz ] can be used to recover a y estimate (and its covariance) 
when an estimate for x (and its covariance) are determined. (See example 
II. 4) . 

Still another application related to the combining of data sets involves 
the combining of SRIF triangular data arrays. One might encounter such prob- 
lems when combining data from different space missions (that involve common 
parameters) or one might choose to process data of each type* or tracking, 
station separately and then combine the resulting SRIF arrays . Triangular 
arrays can be combined using subroutine TTHH, assuming that subroutines 
R2A, THH and R2RA have been used previously to formulate a common parameter 
set for each of the sub problems. 

II. 7 Batch Sequential White Noise 

It is not uncommon to have a problem where each data set contains a 
set of parameters that ^.pply only to that set and not to any other, viz. 
the data is of the form 


A.x+B.y. - z . -v. 
3 3 3 3 3 


3 1, . . . ,N 


where there is generally a priori information on the vector y. variables. 
Rather than form a concatenated state vector composed of x, y^, . . . ,y 
which might create a problem involving exhorbitant amounts of storage and 
computation we solve the problem as follows. Apply subroutine THH to 
[B^:A^:z^3> with the corresponding R initialized with the y^ a priori. The 
resulting SRIF array is of the form 

* 

viz. range, doppler, optical, etc. 
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Copy the top N rows if one will later want an estimate or covariance of 
y l 

the y parameters. Apply subroutine TZERO to zero the top N rows and 
1 y l 

usmg subroutine R2RA set in the y^ a priori . This SRIF array is now 

ready to be combined with the second set of data z^] and the procedure 

repeated. 

A somewhat analogous situation is represented by the class of problems 
that involve noisy model variations, i.e., the state at step j+1 satisfies 


x. = x. + G. 
1+1 3 J 


w. 

3 


where matrix G is defined so that w. is independent of x and w eNCO.O ). 

3 3 3 3 j 

Models of this type are used to reflect that the problem at hand is not 


truly one of parameter estimation, and that some (or all) of the components 
vary in 'a random (or at least unknown) manner that is statistically 
bounded. To solve this problem in a SRIF formulation suppose that a priori 
f° r x and w are written in data equation form (cf ref. [3]), 


R.x. = z. - v. 
3 3 3 3 


v j eN(0,I) 


Q“. 1/2 w. = 0 - 
3 3 3 


vf w) eN(0,I ) 
3 n 

w 


where Q. is a Cholesky factor of Q. that is obtainable from C0V2RI. 
J 3 

these two equations with the one for x gives 

3+1 


In this example it is assumed that all of the yj variables have the same 
dimension. This assumption, though not essential, simplifies our description 
of the procedure. 
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I 

n 

w 

0 


i 

i 


0 


(w) 

v; 

3 

-R.G.Q 3 ? 
3 3 3 

R. 

J 


1 


z . 
3 


v. 

J 


where = w.. . This is the equation to be triangularized with subroutine 

THH, i.e., 


Dim w Dim x 1 


Dim w { 

I 

n 

w 

0 

0 

THH 

3 

r(wx> 

J 

w 
z . 
J 

Dim x { 

k 

-R.G.Q? 
_ J 3 3 

R. 

3 

z 


0 

R. ., 
3+1 

V 1_ 


When the problem is arranged so that Q_. is diagonal one can reduce storage 
and computation. Incidentally, the form of this algorithm allows one to use 
singular Q matrices. 


When the a priori for x^ and Q are given in U-D factored form, 

one can obtain the U-D factors for x. . , as follows: 

3+1 

Let Qj = U (q) D (q) (U^) T (use C0V2UD if necessary) 

Set G = Gj U (q) = [ gl , g n ] , D (q) = Diag(d 1 ,...,d n ) 

w w 


Apply subroutine RANK1 n times, with U n = U. , D n = D. 


w 


o “ u j 5 J 


, RANK1 , 

(u— D) k ; d k , g k (U-D) 


k+1 


i.e. 


-T 


( Wk + d A< - 


k = 1, . . . ,n 


w 


Then U = U , D ... = D 
3+1 n i+l n 
J w J w 


SSSS ® 5 
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Certain filtering problems involve dynamic models of the form 


x. , - = x. + G. w. 
1+1 ii J i 


Given an estimate for x. , x,, the predicted estimate for x..,. 

11 • 1 + 1 - 


denoted 


x j+ i is simply 


x _ = X. 

1+1 1 1 

The U-D factors of the estimate error corresponding to the estimate x^ + ^ 
can he obtained using the weighted Gram-Schmidt triangularization subroutine 
[$. V. : G]; Diag (D jf D (q) ) W G S - K TJ j+1 - D j+1 ) 

Subroutine PH1U can be used to construct Note that this matrix multi- 


1 1 

plication updates the estimate too, because it is placed as an addended column 
to the U matrix. 

When the w and associated x terms correspond to a colored noise model, 

p, t1 = mp,+ w., then it is easier and more efficient to use the colored noise 
1+1 1 1 

update subroutine UDCOL. Note that here too the estimate is updated by the 
subroutine calculation because the estimate is an addended column of U. 

II. 8 Miscellaneous Uses of the Various ESP Subroutines 

In certain parameter analyses we may want to reprocess a set of data 
suppressing different subsets of variables . In this case the original data 
should be left unaltered and subroutine A2A1 usdd to copy A into A^, which 
then can be modified as dictated by the analysis . 

Covariance analysis sometimes are initialized using a covariance 
matrix from a different problem (or a different phase of the same problem) . 

In such cases it may be necessary to permute, delete or insert rows and 
columns into the covariance matrix; and that can be achieved using sub- 
routine C2C. 

If a priori for the problem at hand is given as a covariance matrix 

then one can compute the corresponding SRIF or U-D initialization using 
_ 


In statistical notation that is commonly used, one writes 
*(j+l|j) = x(j|j) 
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subroutines C0V2RI or C0V2UD. Of course, if the covariance is diagonal 
the appropriate R and U-D factors can be obtained more simply. To 
convert a priori given in the form of an information matrix to a corres- 
ponding SRIF matrix one applies subroutine INF2R. To display covariance 
results corresponding to the SRIF or U-D filter one can use subroutines 
UTINV, RI2C0V and UD2C0V. The vector stored covariance results can be 
displayed in a triangular format using subroutine TWOMAT. 

Parameter estimation does not, in the main, involve matrix multipli- 
cation. Certain applications, such as coordinate transformations and time 
propagation are important enough to warrant inclusion in the ESP. For that 
reason we have included RA (to post multiply a square root information 
matrix) and PHIU to premultiply a U-covariance factor) . Certain time propa- 
gation problems involve sparse transition matrices, and for this we have 
included the subroutine SFU. Other special matrix products involving tri- 
angular matrices were not included because we have had no need for other 
products (to date) , and they are generally not lengthy or complicated to 
construct. We illustrate this point by showing how to compute z = Rx where 
R is a triangular vector stored matrix and x is an N vector. 


11=0 


DO 2 1=1, N 


SUM=0. 

@SUM is Double Precision 

11=11+1 

@II=(I,I) 

IK=II 


DO 1 K=I,N 


SUM=SUtH-R (IK) *x(K) 

@IK=(I,K) 

1 IK=IK+K 


2 z(I)=SUM 

@z can overwrite x if desired 
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Note that the II and IK incremental recursions are used to circumvent 
the N(N+1) / 2 calculations of IK=K(K-1) /2+I_. 
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III. SUBROUTINE DIRECTORY SUMMARY 


1. A2A1 - (A to Al) 

Reorders the columns of a rectangular matrix A, storing the 
result in matrix Al. Columns can be deleted and new^ columns added. 
Zero columns are inserted which correspond to new column name entries. 
Matrices A and Al cannot share common storage. 

Example III.l 


a 

B 

C 

B 

F 

G 

C 

H 

1 

5 

9 


5 

0 

0 

9 

0 

2 

6 

10 

A2A1 

6 

0 

0 

10 

0 

3 

7 

11 


7 

0 

0 

11 

0 

-4 

8 

12 


8 

0 

0 

12 

0 


A Al 

The new namelist (BFGCH) contains F, G and H as new columns and deletes 
the column corresponding to name a. 


Example, III. 2 

Suppose one is given an observation data file with regression 

coefficients corresponding to a state vector with components say, 

x, y, z, v , v , v and station location errors. Suppose further, 
x y z 

*f* *f* -f* 

that the vector being estimated has components a^, a^, a y } 

x, y, z, v , v , v , GM and station location errors. A2A1 can be used 
x y z 

to reorder the matrix of regression coefficients to correspond to the 
state being estimated. Zero coefficients are set in place for the 
accelerations and GM which are not present in the original file. 


in track and cross track accelerations 
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2 . 


COMBO - (combine R and A namelists) 


The upper triangular vector stored matrix R has its columns 
permuted and is copied into matrix A. The names associated with R 
are to be combined with a second namelist. 

The namelist for A is arranged so that R names not contained in 
the second list appear first (left most). These are then followed by 
the second list. Names in the second list that do not appear in the 
R namelist have columns of zeros associated with them. 

Example III, 3 

NAM2 list 


a 

B 

C 

D 

C 

B 

E 

a 

F 

D 

~ 1 

2 

4 

7 " 


" 4 

2 

0 

1 

0 

7 “ 

0 

3 

5 

8 


5 

3 

0 

0 

0 

8 

0 

0 

6 

9 


6 

0 

0 

0 

0 

9 

0 

0 

0 

10 


0 

0 

0 

0 

0 

10 


R-Vector stored A-Double subscripted 

A principal application of this subroutine is to the problem of 
combining equation sets containing different variables, and automating 
the process of combining name lists. 

3. COVRHO - (covariance to correlation matrix) 

A vector stored correlation matrix, RHO, is computed from an 
input positive semi-definite vector stored matrix, P. Correlations 
corresponding to zero diagonal covariance elements are zero. To econo- 
mize on storage the output RHO matrix can overwrite the input P matrix. 

The principal function of correlation matrices is to expose strong 
pairwise component correlations ( |rHO(IJ) | .LE . 1, and near unity in magni- 
tude) . It is sometimes erroneously assumed that numerical ill— conditioning 
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of the covariance matrix can. be determined by inspecting the correlation 
matrix entries. While it is true that RHO is better conditioned than is 
the covariance matrix, it is not true that inspection of RHO is sufficient 
to detect numerical ill-conditioning. For example, it is not at all 
obvious that the following correlation matrix has a negative eigenvalue. 


RHO = 


1. 0.50001 

1 . 


0. 50001 
-0.50001 

1 . 


4. C0F2RI - (Covariance to R inverse) 

An input positive semi-definite vector stored matrix P is replaced 

T 

by its upper triangular vector stored Cholesky factor S, P = SS . The name 
RI is used because when the input covariance is positive definite, S = R \ 

5. C0V2UD - (Covariance to U-D factors) 

An input positive semi-definite vector stored matrix P is replaced 

T 

by its upper triangular vector stored U-D factors . P = UDU . 


6. C2C - (C to C) 

Reorders the rows and columns of a square (double subscripted) 
matrix C and stores the result back in C. Rows and columns of zeros 
are added when new column entries are added. 

Example III. 4 



A 

B 

r 



r 

P 

B 

Q 

A 

1 

4 

7 


r 

9 

0 

6 

0 

B 

2 

5 

8 

C2C 

p 

0 

0 

0 

0 

r 

3 

6 

9 


B 

8 

0 

5 

0 






Q 

0 

0 

0 

0 


Names P and Q have been added and name A deleted. An important appli- 
cation of this subroutine is to the rearranging of covariance matrices. 
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7. 


INF2R - (Information matrix to R) 


Replaces a vector stored positive semi-definite information matrix 
A by its lower triangular Cholesky factor R ; A = R R. The upper tri- 
angular matrix R is in the form utilized by the SRIF algorithms. The 
algorithm is designed to handle singular matrices because it is a 
common practice to omit a priori information on parameters that are 
either poorly known or which will be well determined by the data. 

8. HHPOST - (Householder orthogonal triangularization by post 
multiplication) 

The input, double subscripted, rectangular matrix W(M,N) (M.LE.N) 
is triangularized, and overwritten, by post-multiplying it by an implicitly 
defined orthogonal transformation, i.e. 

[ W ]T *■[ 0 \S] 


This subroutine is used, in the main, to retriangularize a mapped covari- 


ance square root and to include in the effects of process noise (i.e. 

1/2 1/2 

W = [§ *P :BQ ]) and to compute consider covariance matrix square 

roots (i.e. W = [P 1/2 _ • Sen*? 1 / 2 ]). 

computed y J 

9 . PERMUT 


Reorders the columns of matrix A, storing the result back in A. 
This routine differs from A2A1 principally in that here the result over- 
writes A. PERMUT is especially useful in applications where storage is 
at a premium or where the problem is of a recursive nature. 

10 . PHIU - (PHI (rectangular) * U (unit upper triangular) ) 


[ PHI ] 



= [ PHIU 3 


The matrices PHI and PHIU are double subscripted, and U is vector sub- 
scripted with implicitly defined unit diagonal elements. It is not 
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necessary to include trailing columns of zeros in the PHI matrix; they 
are accounted for implicitly. To economize on storage the output PHIU 
matrix' can overwrite the input PHI matrix. For problems Involving sparse 
PHI matrices it is more efficient to use the sparse matrix multiplication 
subroutine, SFU. When the last column of U contains the estimate, x, the 
last column of W represents the mapped elements of PHI*x. The principal 
use of this subroutine is the mapping of covariance U factors , where P = UDU , 
and estimates. 

11. RA - (R (triangular) * A(rectangular) ) 



Square root information matrix mapping involves matrix multipli- 
cation of the form indicated in the figure, i.e. with the bottom portion 
of A only implicitly defined as a partial identity matrix. Features of 
this subroutine are that the resulting RA matrix can overwrite the input 
A, and one can compute RA based on a trapezoidal input R matrix (i.e. only 
compute part of R*A) . 

12. RANK1 - (U-D covariance factor rank 1 modification) 

Computes updated U-D factors corresponding to a rank 1 matrix 

modification; i.e., given U-D, a scalar c, and vector v, U and D are 
- — — T T T 

computed so that UDU =UDU +cvv. Both c and v are destroyed during 

the computation, and the resultant (vector stored) U-D array replaces 

the original one. Uses for this routine include (a) adding process 

noise effects to a U-D factored Kalman filter; (b) computing consider 

covariances (cf Section II. 5); (c) computing "actual" covariance 

factors resulting from the use of suboptimal Kalman filter gains; and 

(d) adding measurements to a U-D factored information matrix. 
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13. 


RCOLRD - (colored noise inclusion into the SRIF) 


Includes colored noise time updating into the square root infor- 
mation matrix. It is assumed that the deterministic portion of the time 
update has been completed, and that only the colored noise effects are 
being incorporated by this subroutine. The algorithm used is Bierman’s 
colored noise one-component-at-a-time update, cf ref. [3], and updates the 
SRIF array corresponding to the model 


V 


'I 

0 

0 " 


V 


" o' 

p 

= 

0 

M 

0 


p' 

+ 

w j 

X 2 

3+1 

0 

0 

I 


_ X 2_ 

j 

0 


M is diagonal and w^ eH(0,Q). Auxiliary quantities, useful for fixed interval 
smoothing, are also generated. 

14. RINCON - (R inverse with condition number bound, CNB) 

Computes the inverse of an upper triangular vector stored matrix R 
using back substitution. To economize on storage the output result can 
overwrite the input matrix. A Frobenius bound (CNB) for the condition 
number of R is computed too. This bound acts as both an upper and a 
lower bound, because CNB/N £ condition number < CNB. When this bound is 
within several orders of magnitude of the machine accuracy the computed 
inverse is not to be trusted, (viz if CNB ^ 10^ on an 18 decimal digit 
machine R is ill-conditioned) . 

15. RI2C0V - (RI to covariance) 

This subroutine computes sigmas (standard deviations) and/or the 
covariance of a vector stored upper triangular square root covariance 
matrix, RINV (SRIF inverse). The result, stored in C0V0UT (covariance 
output) is also vector stored. To economize on storage, COVOUT can over- 


write RINV. 



16. 


R2A - (R to A) 


The columns of a vector stored upper triangular matrix R are per- 
muted and variables are added and/or deleted. The result is stored in 
the double subscripted matrix A. In other respects the subroutine is 
like A2A1. 

Example III. 5 


a 

B 

C 

D 

E 

E 

F 

C 

B 

" 2 

4 

8 

14 

22" 


22 

0 

8 

4" 

0 

6 

10 

16 

24 


24 

0 

10 

6 

0 

0 

12 

18 

26 

R2A 

26 

0 

12 

0 

0 

0 

0 

20 

28 


28 

0 

0 

0 

0 

0 

0 

0 

30 


30 

0 

0 

0 


R A 

R is vector stored as R = (2,4,6,8,10,12,14,16,18,20,22,24,26,28,30) 
with namelist (a,B,C,D,E) associated with it. Names a and D are 
not included in matrix A, and a column of zeros corresponding to name 
F is added. 

One trivial, but perhaps useful, application is to convert a 
vector stored matrix to a double subscripted form.^ R2A is used, most 
often when one wants to rearrange the columns of a SRIF array so that 
reduced order estimates, sensitivities, etc. can be obtained; or so that 
data sets containing different parameters can be combined. 


t 


see also the aside in the introduction 
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17. R2RA - (Triangular block of R to triangular block of RA) 

A triangular portion of the vector stored upper triangular matrix R 
is put into a triangular portion of the vector stored matrix RA. The 
names corresponding to the relocated block are also moved. R can coin- 
cide with RA. 



Note that an upper left triangular submatrix can slide to any lower 
position along ,the diagonal, but that a submatrix moving up must go 
the upper leftmost corner. Upper shifting is used when one is 
interested in that subsystem; and the lower shifting is used, for 

example, when inserting a priori information for consider analyses. 
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18. RUDR - (SRIF R converted to U-D form or vice versa) 

A vector stored SRIF array is replaced by a vector stored U-D 
form or conversely. A point to be noted is that when data is involved 
the right side of the SRIF data equation transforms to the estimate in 
the U-D array. 

19. SFU - (Sparse F*U(Unit upper triangular)) 



A sparse F matrix, with only its nonzero elements recorded, multiplies 
U which is vector stored with implicit unit diagonal entries. When the 
input F is sparse this routine is very efficient in terms of storage and 
computation. When the last column of U contains the estimate, x, the last 
column of FU represents elements of the mapped estimate F * x. 

20. TDHHT - (Two dimensional Householder Triangularization) 

Implicitly defined Householder orthogonal transformations are used 
to triangularize an input two dimensional rectangular array, S(M,N)* 
Computation can be reduced if S starts partially triangular; ' 

S = 

JSTART 

Further, the algorithm implementation is such that (a) maximum trian- 
gularization is achievable 

when M.LT.N S -> 
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when M. GT . N S 


and finally when an intermediate form is desired 



V 


s 

0 

^ - 





JSTOP 

This subroutine can be used to compress overdetermined linear systems of 
equations to triangular form (for use in least squares analyses) . The 
chief application, that we have in mind, of this subroutine, is to the 
matrix triangularization of a "mapped" square root information matrix. 

This subroutine overlaps to a large extent the subroutine THH which 
utilizes vector stored, single subscripted, matrices. This latter rou- 
tine when applicable is more efficient. The triangularization is adapted 
from ref. [1], 

21. THH - (Triangular Householder data packing) 

An upper triangular vector stored matrix R is combined with a 
rectangular doubly subscripted matrix A by means of Householder orthogonal 
transformations. The result overwrites R, and A is destroyed in the process. 
This subroutine is a key component of the square root information sequential 
filter, cf ref. [3]. 


THH 


* 

The elements are not explicitly set to zero. 





32 



22. 


TTHH - (Two triangular arrays are combined using Householder 
orthogonal transformations) 

This subroutine combines two single subscripted upper triangular 
SRIF arrays, R and RA using Householder orthogonal transformations. The 
result overwrites R. 



23. TWOMAT - (Two dimensional print of a triangular matrix) 

Prints a vector stored upper triangular matrix, using a matrix 

format. 


Example III, 7 

R(10) = (2,4,6,8,10,12,14,16,18,20) with associated namelist 
(A,B,C,D) is printed as 

A B C D 

A 2 4 8 14 

B 6 10 16 

C 12 18 

D 20 

(The numbers are printed as 7 columns of 8 significant 
floating point digits or 12 columns of 5 significant floating 
point digits.) 

To appreciate the importance of this subroutine compare the vector 
R(10) with the double subscript representation. 

The elements are not explicitly set to zero. 
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24. TZERO - (Zero a horizontal segment of a vector stored upper 


triangular matrix) 

Upper triangular vector stored matrix R has its rows between ISTART 
and IFINAL set to zero. 

Example III. 8 

To zero rows 2 and 3 of R(15) of example III. 5 

R(15) = (2,4,6,8,10,12,14,16,18,20,22,24,26,28,30) is transformed to 
R(15) = (2,4,0,8,0,0,14,0,0,20,22,0,0,28,30) i.e.. 


2 

4 

8 

14 

22 


2 

4 

8 

14 

22 

0 

6 

10 

16 

24 


0 

0 

0 

0 

0 

0 

0 

12 

18 

26 

TZERO. 

0 

0 

0 

0 

0 

0 

0 

0 

20 

28 


0 

0 

0 

20 

28 

0 

0 

0 

0 

30 


0 

0 

0 

0 

30 


R-vector stored R-vector stored 


25. UDCOL - (U-D covariance factor colored noise update) 

This subroutine updates the U-D covariance factors corresponding 
to the model 



0 0 


i 

X 

1-1 i 


0 

M 0 


P 

+ 

w. 





1 

0 I 


x„ 


0 



2 



• 


- — 

j 



where M is diagonal and w^, £N(0,Q). The special structure of the transi- 
tion and process noise covariance matrices is exploited, cf Bierman, [3]. 
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26. UDMEAS — (TJ— D Measurement Update) 

Given the U-D factors of the a priori estimate error covariance 
and the measurement, z = Ax + V this routine computes the updated estimate 
and U-D covariance factors, the predicted residual, the predicted residual 
variance, and the normalized Kalman gain. This is Bierman’s U-D measure- 
ment update algorithm, cf [3]. 

27. UP 2 CQV - (U-D factors to covariance) 

The input vector stored U-D matrix (diagonal D elements are stored 

as the diagonal entries of U) is replaced by the covariance P, also vector 
T 

stored, P = UDU . P can overwrite U to economize on storage. 

28. UD2SIG - (U-D factors to sigmas) 

Standard deviations corresponding to the diagonal elements of the 
covariance are computed from the U-D factors. This subroutine, a restricted 
version of UD2C0V can print out the resulting sigmas and a title. The 
input U-D matrix is unaltered. 

29. UTINV - (Upper triangular matrix inversion) 

An upper triangular vector stored matrix KIN (R in) is inverted 
and the result, vector stored, is put in ROUT (R out). ROUT can overwrite 
RIN to economize on storage. If a right hand side is included and the 
bottommost tip of RIN has a -1 set in then ROUT will have the solution in 
the place of the right hand side. 
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30. UTIROW - (Upper triangular inversion, inverting only the upper rows) 


INPUT OUTPUT 


R 

X 

R 

xy 

UTIROW 

R 1 

X 

-1 -d -1 

” R x R xy y 

0 

R" 1 

0 

r" 1 


y 



y 


An input vector stored R matrix with its lower left triangle assumed to 
have been already inverted is used to construct the upper rows of the 
matrix inverse of the result. The result, vector stored, can overwrite 
the input to economize on storage. 

If the columns comprising R^ represent consider terms then taking 
R ^ as the identity gives the sensitivity on the upper right portion of 
the result. If R ^ = Diag (0 ,...,0 ) then the upper right portion of 

y y n y 

the result represents the perturbation. Note that if z (the right hand 
side of the data equation) is included in R^ then taking the corres- 
ponding R ^ diagonal as -1 results in the filter estimate appearing 
as the corresponding column of the output array. When n^ is zero this 
subroutine is algebraically equivalent to UTINV. The subroutines differ 
when a zero diagonal is encountered. UTINV gives the correct inverse 
for the columns to the left of the zero element, whereas UTIROW gives 
the correct inverse for the rows below the zero element. 
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ORIGINAL PAGE IS 
OP POOR QUALITY 


31. WGS - (Weighted Gram-Schmidt U-D matrix triangularization) 

An input rectangular (possibly square) matrix W and a diagonal 

weight matrix, D , are transformed to (U-D) form; i.e,, 
w 

T T 

S D W = UDU 
w 

where U is unit upper triangular and D is diagonal. The weights D are 

w 

assumed nonnegative, and this characteristic is inherited by the 
resulting D. 
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IV . SUBROUTINE DIRECTORY USER DES CRIPTION 


1. A2A1 (A to Al) 


Purpose 

To rearrange the columns of a namelist indexed matrix to 


conform to a desired namelist. 


CALL A2A1 (A > IA,IR > LA > NAMA > A1 > IA1 > LA1 > NAMA1) 


Argument Definitions 


A(IR,LA) 


Input rectangular matrix 


IA 


Row dimension of A, IA.GE.IR 


IR 


Number of rows of A that are to be 
arranged 


LA Number of columns in A; this also 

represents the number of parameter 
names associated with A 


NAMA(LA) 


Parameter names associated with A 


A1(IR,LA1) Output rectangular matrix 

IA1 Row dimension of Al, IA1.GE.IR 


LAI Number of columns in Al ; this also 

represents the number of “parameter 
names associated with Al 


NAMAl(LAl) 


Input list of parameter names to be 
associated with the output matrix Al 


Remarks and Restrictions 


Al cannot overwrite A. This subroutine can be used to add 
on columns corresponding to new names and/or to delete variables 


from an array. 

Functional Description 

The columns of A are copied into Al in an order corresponding 

to the NAMA1 parameter namelist. Columns of zeros are inserted 

♦ 

in those Al columns which do not correspond to names in the input 
parameter namelist NAMA. 
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2. COMBO (Combine parameter namelists) 

Purpose 

To rearrange a vector stored triangular matrix and store 
the result in matrix A. The difference between this subroutine 




OF 


and R2A is that there the namelist for A is input; here it is 
determined by combining the list for R with a list of desired names. 

CALL COMBO (R,L1,NAM1,L2,NAM2,A,IA,LA,NAMA) 

Argument Definitions 


R(L1*(L1+1) /2) Input vector stored upper triangular matrix 

LI No. of parameters in R (and in NAMl) 

KAMI (LI) Names associated with R 


L2 


No . of parameters in NAM2 


NAM2 (L2) 


A (LI, LA) 


Parameter names that are to be combined 
with R (NAMl list); these names may or 
may not be in NAMl 

Output array containing the rearranged 
R matrix Ll.LE.IA 


IA 


Row dimension of A 


LA No. of parameter names in NAMA, and the 

column dimension of A. LA=Ll + L2 - 
No. names common to NAMl and NAM2; LA 
is computed and output 

NAMA(LA) Parameter names associated with the out- 

put A matrix ; consists of names in NAMl 
which are not in NAM2, followed by NAM2 


Remarks and Restrictions 


The column dimension of A is a result of this subroutine. 

To avoid having A overwrite neighboring arrays one can bound the 
column dimension of A by L1 + L2. 
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Functional Description 

First the KAMI and NAM2 lists are compared and the names 
appearing in NAM1 only have their corresponding R column entries 
stored in A (e.g. if NAML(2) and NAM1(6) are the only names not 
appearing in the NAM2 list then columns 2 and 6 of R are copied 
into columns 1 and 2 of A) . The remaining columns of A are 
labeled with NAM2. The A namelist is recorded in NAMA. The 
NAM1 list is compared with NAM2 and matching names have their R 
column entries copied into the appropriate columns of A. NAM2 
entries not appearing in NAMl have columns of zero placed in A. 
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3. COVRHO (Covariance to correlation matrix, RHO) 

Purpose 

To compute the correlation matrix RHO from an input covariance 

matrix COV. Both matrices are upper triangular, vector stored and 
the output can overwrite the input. 

| CALL COVRHO (COV, N, RHO, V) 

Argument Definitions 


COV (N* (N+l) /2) 
N 

RHO (N* (N+l)/2) 
V (N) 


Input vector stored positive semi-definite 
covariance matrix 

Model dimension, N.GE.l 

Output vector stored correlation matrix 

Work vector 


Remarks 

No test for non-negativity of the input matrix is made. 


Correlations corresponding to negative or zero diagonal entries 


are set to zero, as is the diagonal output entry. 

Functional Description 

V(I). = 1//C0V(I,I) if COV(I,I),GT.O and 0. otherwise 
RH0(I,J) = C0V(I,J)*V(I)*V(J) 

The subroutine employs, however, vector stored COV and RHO matrices. 
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4. C0V2RI (Covariance to Cholesky Square Root, RI) 

Purpose 

To construct the upper triangular Cholesky factor of a positive 
semi-definite matrix. Both the input covariance and the output 
Cholesky factor (square root) are vector stored. The output 
overwrites the input. Covariance (input) = (CF)*(CF)**T 
(output CF = Rinverse). If the input covariance is singular, the 
output factor has zero columns. 

CALL C0V2RI(CF,N) 

Argument Definitions 

CF(N* (N+l)/2) Contains the input vector stored 

covariance matrix (assumed positive 
definite) and on output it contains 
the upper triangular Cholesky factor 

N Dimension of the matrices involved, N.GE.2 

Remarks and Restrictions 

No check is made that the input matrix is positive semi-definite. 
Singular factors (with zero columns) are obtained if the input is 
(a) in fact singular, (b) ill-conditioned, or (c) in fact indefinite; 
and the latter two situations are cause for alarm. Case (c) and 
possibly (b) can be identified by using RI2C0V to reconstruct the 
input matrix. 

Functional Description 

An upper triangular Cholesky reduction of the input matrix is 

implemented using a geometric algorithm described in Ref. [3]. 

CF(input) = CF (output) *CF (output) ^ 

At each step of the reduction diagonal testing is used and negative 

terms are set to zero. 
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5. COV2UD (Covariance to TID factors) 

Purpose 

To obtain the U-L factors of a positive semi-definite matrix. 

The input vector stored matrix is overwritten by the output U-D 
factors which are also vector stored. 

CALL C0V2UD(U,¥j 

Argument Definitions 

U(N*(N+l)/2) Contains the input vector stored covari- 

ance matrix ; on output it contains the 
vector stored U-D covariance factors. 

N Matrix dimension, N,GE,2 

Remarks and Restrictions 

No checks are made in this routine to test that the input U matrix 
is positive semi-definite. Singular results (with zero columns) are 
obtained if the input is (a) in fact singular, (b) ill-conditioned, 
or (c) in fact indefinite; and the latter two situations are cause for 
alarm. Case (c) and possibly case (b) can be identified by using UD2-> 
COV to reconstruct the input matrix. Note that although indefinite 
matrices have U-D factorizations, the algorithm here applies only to 
matrices with non-negative eigenvalues. 

Functional Description 

An upper triangular U-D Cholesky factorization of the input matrix 

is implemented using a geometric algorithm described in Ref. [3]. 

T 

U (input) = U*D*U , U-D overwrites the input U 

at each step of the reduction diagonal testing is used to zero negative 
terms. 
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6. C2C (C to C) 

Purpose 

To rearrange the rows and columns of C, from NAMl order to NAM2 
order. Zero rows and columns are associated with output defined names 
that are not contained in NAML. 


CALL C2C(C,IC,L1,NAM1,L2,NAM2) 

Argument Definitions 

C (LI ,Ll) Input matrix 

IC Row dimension of C 

IC.GE.L = MAX(Ll,L2) 


LI No . of parameter names associated with 

the input C 

NAKL(L) Parameter names associated with C on input. 

(Only the first LI entries apply to the 
input C) 

L2 No. of parameter names associated with the 

output C 


NAM2(L2) 


Parameter names associated with the output C 


Remarks and Restrictions 

The NAM2 list need not contain all the original NAMl names and 
LI can be .GE. or .LE. L2. The NAMl list is used for scratch and 
appears permuted on output. If L2.GT.L1 the user must be sure that 
NAMl has L2 entries available for scratch purposes . 

Functional Description 

The rows and columns of C and NAMl are permuted pairwise to get 
the names common to NAML and NAM2 to coalesce. Then the remaining rows 
and columns of C(L2,L2) are set to zero. 
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7 . HHPOST (Householder Post Multiplication Triangularization) 

Purpose 

To employ Householder orthogonal transformations to triangularize 
an input rectangular W matrix by post multiplication, i.e. 



This algorithm is employed in various covariance square root updates . 

CALL HHPOST (S,W,MROW, NROW, NCOL, V) | 

Argument Definitions 
S (NROW* (NROW+1) /2) 

W (NROW, NCOL) 

MROW 
NROW 

NCOL 
V(NCOL) 

Functional Description 

Elementary Householder transformations are applied to the rows of W 
in much the same way as they are applied to obtain subroutine THH. The 
orthogonolization process is discussed at length in the books by Lawson 
and Hanson [1] and Bierman [3]. 


Output upper triangular vector stored 
square root matrix 

Input rectangular square root covariance 
matrix (W is destroyed by computations) 

Maximum row dimension of F 

Number of rows of W to be triangular! zed 
and the dimension of S (NRQW.GE.2) 

Number of column of W (NCOL. GE. NROW) 

Work vector 
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8. INF2R (Information matrix to R) 

Purpose 

To compute a lower triangular Cholesky factorization of an 
input positive semi-definite matrix. The result transposed, is 
vector stored; this is the form of an upper triangular SRIF matrix. 

CALL INF2R(R,N) 

Argument Definitions 

R (N* (N+l) / 2) Input vector stored positive semi- 

definite (information) matrix; on output 
it represents the transposed lower 
triangular Cholesky factor (i.e. the SRIF 
R matrix) 

N Matrix dimension, N.GE.2 

Remarks and Restrictions 

No checks are made on the input matrix to guard against negative 
eigenvalues of the input, or to detect ill-conditioning. Singular 
output matrices have one or more rows of zeros . 

Functional Description 

A Cholesky type lower triangular factorization of the input matrix 

is implemented using the geometric formulation described in Ref. [3]. 

T 

R(input) = [R (output) ] * [R(output) ] 

At each step of the factorization diagonal testing is used to zero columns 
corresponding to negative entries. The result is vector stored in the 
form of a square root information matrix as it would be used for SRIF 
analyses . 
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9 . PERMUT (Permute A) 

Purpose 

To rearrange the columns of a namelist indexed matrix to conform 
to a desired namelist. The resulting matrix is to overwrite the input. 


CALL PERMUT (A, IA, IR,L1 ,NAM1 ,L2 ,NAM2^ 


Argument Definitions 
A(IR,L) 

IA 

IR 

LI 

NAM1(L) 

L2 

NAM2 


Input rectangular matrix, L = max(Ll,L2) 

Row dimension of A, IA.GE.IR 

Number of rows of A that are to be 
rearranged 

Number of parameter names associated with 
the input A matrix 

Parameter names associated with A on input 
(only the first LI entries apply to the 
input A) 

Number of parameter names associated with 
the output A matrix 

Parameter names associated with the output A 


Remarks and Restrictions 


This subroutine is similar to A2A1; but because the output matrix 
in this case overwrites the input there are several differences. The 
NAMl vector is used for scratch, and on output it contains a permuta- 
tion of the input NAMl list. The user must allocate L = max(Ll,L2) 
elements of storage to NAMl. The extra entries, when L2>L1, are 
used for scratch. 

Functional Description 

The columns of A are rearranged, a pair at a time, to match the 
NAM2 parameter namelist. The NAMl entries are permuted along with the 
colu mn s, and this is why dim (NAMl) must be larger than LI (when L2>L1) . 
Columns of zeroes are inserted in A which correspond to output names 
that do not appear in NAMl. 
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10. PHIU (PHI-rectangular*U-unit upper triangular) 


Purpose 

To multiply a rectangular two dimensional matrix PHI by a unit 
upper triangular vector stored matrix U, and store the result in 
PHIU. The PHIU matrix can overwrite PHI to economize on storage. 


[PHI] 



= [PHIU] 


CALL PHIU (PHI , MAXPHI , IRPHI , JCPHI ,U ,N ,PHIU,MPHIU) 


Argument Definitions 

PHI (IRPHI, JCPHI) 

MAXPHI 

IRPHI 

JCPHI 

U (N* (N+l) / 2) 

N 

PHIU ( IRPHI, N) 
MPHIU 


Input rectangular matrix IRPHI. LE MAXPHI 

Row dimension of PHI 

number of rows of PHI 
number of columns of PHI 

unit upper triangular vector stored matrix 
U-matrix dimenstion, JCPHI. LE.N 
output result PHI*U , PHIU can overwrite PHI 
row dimension of PHIU 


Remarks and Restrictions 


If JCPHI. LT.N it is assumed that there are implicitly defined 
trailing columns of zeros in PHI. The unit diagonal entries of U 
are implicit, i.e. the diagonal U entries are not explicitly used. 
Functional Description 
PHIU = PHI*U 
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11. RA (R-upper triangular*A-rectangular) 

Purpose 

To post multiply a vector stored triangular matrix., R, by a 
rectangular matrix A, and if desired to store the result in A. 

[ a ] - [ “ ] 


CALL RA(R,N,A,MAXA, IA, JA,RA,MAXRA, IRA) 


Argument Definitions 

upper triangular, vector stored input 
order of R 

Input rectangular right multiplier matrix 

Row dimension of input A matrix 

Number of rows of A that are input 

Number of columns of A 

Output resulting rectangular matrix 
RA can overwrite A 

Row dimension of RA 

Number of rows in the output result 
(IRA.LE .MAXRA) 

Functional Description 

The first IRA rows of the product R*A are computed using the 
vector stored input matrix R, and the output can, if desired, 
overwrite the input A matrix. When N.GT.IA (i.e. there are more 
columns of R than rows of A) then it is assumed that the bottom 
N-IA rows of A are implicitly defined as a partial identity matrix, i.e. 


R(N*(N+l)/2) 

N 

A(IA,JA) 

MAXA 

IA 

JA 

RA(IRA, JA) 

MAXRA 

IRA 


(Input) 
0 I 


}lA 

}N-IA 
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12. RANK1 (Stable U-D rank one update) 


Purpose 

T T 

To compute the (updated) U-D factors of UDU 4- GW . 

CALL RANK1 (UIN ,UOUT ,N , C , V) 

Argument Definitions 
UIN (N* (N+l) / 2) 

UOUT (N* (N+l) /2) 

N 
C 

V (N) 

Remarks and Restrictions 

If C negative is used the algorithm is numerically unstable, 
and the result may be numerically unreliable. Singular U matrices 
are allowed, and these can result in singular output U Matrices. 
The code switches from a 1-multiply to a 2-multiply mode at a key 
place, based upon a 1/16 comparison of input to output D values. 
Also, there is provision made to supply a machine accuracy epsilon 
when single precision is specified. 

Functional Description 

This rank one modification is based on a result published by 
Agee and Turner (1972) , White Sands Missile Range Tech. Report 
No. 38 and improved on using a numerical stabilization idea due 
to Gentlemen (1973). The algorithm is derived in the chapter. 


Input vector stored positive semi- 
definite U-D array (with the D entries 
stored on the diagonal of U) 

Output vector stored positive (possibly) 
semi-definite U-D result, UOUT=UIN is 
allowed . 

Matrix dimension, N.GE.2 

Input scalar, which should be non-negative. 
C is destroyed by the algorithm. 

Input vector for the rank one modification. 
V is destroyed by the algorithm. 
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T 

"UDU Covariance Factorization For Kalman Filtering," C. L. Thornton, 
G. J. Bierman, Vol. XVI of Advances in Control of Dynamic Systems, 
Academic Press, to appear 1979. 
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13. RCOLRD (Colored noise time update of the SRIF R matrix) 

Purpose 

To include colored noise time updating into the square root 
information matrix. It is assumed that the deterministic portion 
of the time update has been completed, and that only the colored 
noise effects are being incorporated by this subroutine. 


CALL RCOLRD ( S ,MAXS , IRS , JCS , NP STRT , NP , EM , RW , ZW , V , SGSTAR) 


Argument Definitions 


S (IRS, JCS) 


MAXS 

IRS 

JCS 


NPSTRT 


Input rectangular portion of the square 
root information matrix corresponding to 
the nonconstant paramters. It is assumed 
that estimates are included, i.e. the last 
column represents the "right hand side",Z, 

(but see JCS description) . S also houses the 
time updated array, and if there is smoothing 
there are UP extra rows adjoined to S. 

Row dimension of S. If smoothing calculations 
are to be included then MAXS.GE.IRS+NP. 

The number of rows of S, i.e. the number of 
nonconstant parameters (including colored 
noise variables). IRS.GE.2 

The number of columns of S. If the vector 
ZW is zero, then the right hand side of 
transformed estimates need not be included. 

Location of the first colored process noise 
variable. 


NP 


The number of colored noise variables 
contiguous to and following the first. 


EM(NP) 


Vector of exponential colored noise multipliers 
(EM = exp (-DT/TAU)) 


RW (NP) 


Vector of positive reciprocal colored process 
noise standard deviations, i.e. 


p . , = exp (-DT/t ) * p . + w . 

*J+1 F 33 


Rw = 1/0 


w 
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ZW(NP) Vector of normalized process noise a priori 

estimates. ZW is generally zero. 

V(IRS) Work vector. 

SGSTAR(NP) Vector of smoothing coefficients. Needed 

only if smoothing is to be done. 

Remarks and Restrictions 

There are three lines of code associated with smoothing, and 
these are commented out of the nominal case. Therefore, if smoothing 
is contemplated the comments must be removed. The vector SGSTAR is 
involved only with smoothing. Last note: for smoothing, be sure 

that S has NP extra rows to house the smoothing coefficients. 

The ZW vector is generally zero. If ZW = 0 one has the option 
of doing covariance only analyses and the last column of S (the 
right hand side of normalized estimates) can be omitted. 

Because of the large number of arguments appearing in this 
subroutine, and because almost all of them are constant (i.e. with 
succeeding calls only S, and possible EM, RW, ZW and SGSTAR change) 
for a given problem, it is suggested that one a) introduce COMMON, 
b) use this as an internal subroutine, or c) write in-line code. 
Functional Description 

The model is 


X, 


1 

H 

O 

O 




0 

1 




1 



p . 

= 

0 M 0 


P 

4* 

w. 

3 

X 2 


0 0 1 


X 2 


0 



— — 




__ 


1_ L_ i_ — i 


}NPSTRT-1 

}NP 

}N- (NPSTRT-l+NP) 


where M is diagonal, with NP non-negative entries and w_. is a white 

- -1 -T 

noise process with w. €N(w, Q) , Q = R R . The algorithm is based 

2 w w 

on Bierman ’ s one component-at-a-time SRXF time update which economizes 



on storage and computation (see Bierman-Factor-ization Methods for 
Discrete Sequential Estimation, Academic Press 1977). 

When smoothing is contemplated, there is output a vector 0*(NP) 
and a matrix S*(NP,N+1): S* occupies the bottom NP rows of the 
output S matrix. Smoothed estimates of the p terms can be obtained 
from the a* and S* terms as follows : 

Let X* be the previously computed estimates of the N filter 

parameters, then for d[ = NP, NP-l,...l recursively compute 

N 

X*(NSTRT + J-l) := (S*(J, N+l) - £ S*(J,K)X*(K))/a*(j) 

K=1 

Note that the symbol " :==" means is replaced by, so that the old 
values of X*, on the right side, are over-written by the new 
smoothed colored noise estimates. Smoothed covariances can be 
obtained from the S* and o* terms as well, but we do not go into 
detail here; the reader is directed to chapter 10 of the Bierman 
reference . 
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14. RINCON (R inverse with condition number bound) 

Purpose 

To compute the inverse of an upper triangular vector stored 
triangular matrix, and an estimate of its condition number. 


CALL RINCON (RIN, N, ROUT, CNB) 


Argument Definitions 
RIN(N*(N+l)/2) 

N 

ROUT (N*(N+l)/2) 
CNB 


Input vector stored upper triangular matrix 

Matrix dimension, N.GE.2 

Output vector stored matrix inverse 
(RIN *= ROUT is permitted) 

Condition number bound. If k is the 
condition number of RIN, then 
CNB/N.LE.k.LE CNB 


Remarks and Restrictions 

The condition number bound, CNB serves as an estimate of the actual 
condition number. When it is large the problem is ill-conditioned. 
Functional Description 

The matrix inversion is carried out using a triangular back 
substitution. If any diagonal element of the input R matrix is 
zero the condition number computation is aborted. When the first 
zero occurs at diagonal k the matrix inversion is carried out only 
on the first k-1 columns. The condition number bound is computed 
as follows: 

NTOT 

F.NORM R = E R(J) 2 

J=1 

NTOT 


F.NORM R 


_1 - X) 


J=1 
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where NTOT = N*(N+l)/2 is the number of elements in the vector stored 
triangular matrix. The condition number bound, CNB, is given by 
CNB = (F . NORM R * F.NORM R -1 ) 1 / 2 
F.NORM is the Erobenius norm, squared. The inequality 

CNB/N ^ condition number R < CNB 

is a simple consequence of the Frobenius norm inequalities given in 
Lawson-Hanson "Solving Least Squares," page 234. 
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15 . RI2COV (RI Triangular to covariance) 

Purpose 

To compute the standard deviations, and if desired, the 
covariance matrix of a vector stored upper triangular square root 
covariance matrix. The output covariance matrix, also vector 
stored, can overwrite the input. 


CALL RI 2C0V (RINV , N , S IG , COVOUT , KROW , KCOL) 


Argument Definitions 

RINV (N* (N+l) /2 Input vector stored upper triangular 

covariance square root (RINV=Rinverse 
is the inverse of the SRIF matrix) . 

N Dimension of the RINV matrix 


SIG(N) Output vector of standard deviations 

COVOUT (N* (N+l) /2) Output vector stored covariance matrix 

(COVOUT = RINV is allowed) 


f .GT.O 


KROW 


1 .LT.O 


(_ .EQ.O 


Computes the covariance and sigmas 
corresponding to the first KROW variables 
of the RINV matrix 

Computes only the sigmas of the first 
(KROW) variables of the RINV matrix. 

No covariance, but all sigmas (e.g. use 
all N rows of RINV) 


KCOL 


Number of columns of COVOUT that are 
computed. If KCOL.LE.O, then KCOL = KROW. 


Remarks and Restrictions 


Replacing N by |KROW| corresponds to computing the covariance 
of a lower dimensional system. 

Functional Description 
COVOUT=RINV*RINV**T 
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16. R2A (R to A) 

Purpose 

To place the upper triangular vector stored matrix R into the 
matrix A and to arrange the columns to match the desired NAMA para- 
meter list. Names in the NAMA list that do not correspond to any 
name in NAMR have zero entries in the corresponding A columns. 


CALL R2A(R,LR,NAMR,A,IA,LA,NAMA) 


Argument Definitions 
R (LR* (LR+1) / 2 ) 
LR 

NAMR(LR) 

A(LR,LA) 

IA 

LA 

NAMA (LA) 


Input upper triangular vector stored array 

No. of parameters associated with R 

Parameter names associated with R 

Matrix to house the rearranged R matrix 

Row dimension of A, IA.GE.LR. 

No, of parameter names associated with the 
output A matrix. 

Parameter names for the output A matrix. 


Functional Description 

The matrix A is set to zero and then the columns of R are copied 
into A. 
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17 . R2RA (Permute a subportion R of a vector stored triangular matrix) 

A 

Purpose 

To copy the upper left (lower right) portion of a vector stored 
upper triangular matrix R into the lower right (upper left) portion of 
a vector stored triangular matrix RA, 

CALL R2RA(R,NR,NAM,RA,NRA,NAMA) 


Argument Definitions 
R(NR*(NR+l)/2) 
NR 

NAM (NR) 


Input vector stored upper triangular matrix 

Dimension. of vector stored R matrix"^ 

Names associated with R. 


RA(NRA*(NRA+1) /2) Output vector stored upper triangular matrix 

NRA If NRA= 0 on input, then NAMA(l) should have 

the first name of the output namelist. In 
this case the number of names in NAMA, NRA, 
will be computed. The lower right block of 
R will be the upper left block of RA. 

If NRA = last name of the upper left block 
that is to be moved then this upper block 
is to be moved to the lower right corner 
of RA. When used in this mode NRA=NR on 
output’! 


NAMA (NRA) 


Names associated with RA. Note that NRA 
used here denotes the output value of NRA. 


Remarks and Restrictions 

RA and NAMA can overwrite R and NAM. The meaning of the NRA= 0 

option is clarified by the following example: 

C D E 

— » 

12 18 26 
20 28 
30 

_ 

RA 

’’'see the concluding paragraph of Remarks and Restrictions 


A 

B 

C 

D 

E 

2 

4 

8 

14 

22* 


6 

10 

16 

24 



1 12 
I 
1 

18 

26 



20 

28 

- 


1 

R 


30 


INPU T 
NR - 5 

NAM = ’AVB’. 'CVDVE' 
NRA = 0 
NAMA(l) « ’C’ 

R 


OUTPUT 

NAMA = ’C\ *D’ , ’E» 
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When -NBA = 0 and NAMA(l) = 'C* we are asking that the lower triangular 
portion of R, beginning at the column labeled C, be moved to form the 
first (in this case 3) columns of RA, Incidently, RA could have 
additional columns; these columns and their names would be unaltered 
by the subroutine. 

The meaning of the other NRA option is illustrated by the following 
example; 


A B C t D E 

i 1 

A B I A B C h 

r* 1 



2 4 8 1 14 22 


2 4 8 14 22 

6 10 | 16 24 


6 10 16 24 

1 

12 | 18 26 

► 

[2 4 8 

20 28 


1 6 10 

30 


1 12 


R R 


INPUT 
NR = 5 

NAM = ’A\ f B VCVDVE* 
NRA = T C’ 

R 

OUTPUT 
NRA = 5 

NAMA(3-5) = 'A* , ’B' , ! C ' 
RA 


When NRA = T C* we are asking that the upper left block of R, up to the 
column labeled C, be moved to the lower right poriton of RA and the cor- 
responding names be moved too. If RA overwrites R, as in the example, 
then the first two rows of R remain unchanged and since NAMA overwrites 
NAM, the labels of the first two columns remain unaltered. 

The remark that NRA=NR on output means, in this example, that the 
column with name C in R is moved over to column 5. If one wanted to 
slide the upper left triangle corresponding to names ABC of R to columns 
7-9 of an RA matrix (of unspecified dimension, >9), then one should set 
NR=9 in the subroutine call. Thus NR, when used in this sliding down 
the diagonal mode, does not represent the dimension of R; but indicates 
how far the slide will be. 
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18. RUDR (R to U— D or U-D to R) 

Purpose 

To transform an upper triangular vector stored SRIF array to U-D 
form or vice versa. 

CALL RUDR(RIN , N , ROUT , IS ) 


Argument Definitions 

RIN (NBAR* (NBAR+l)/2) Input upper triangular vector stored SRIF 

or U-D array; NBAR = ABS(N) + 1 

ROUT (NBAR* (NBAR+l)/2) Output upper triangular vector stored 

U-D or SRIF array (RIN - ROUT is 
permitted) 


N Matrix dimension, N.GT.O represents an 

R to U-D conversion and N.LT.O represents 
a U-D to R conversion. ABS(N).GE.2 

IS If IS = 0 the input array is assumed not 

to contain a right side (or an estimate) , 
and IS = 1 means an appropriate additional 
column is included. In -the IS = 0 case 
the last column of RIN is ignored and 
NBAR = ABS(N) is used. 


Subroutine used: RINCON 


Functional Description 

Consider the N>0 case. RIN = R is transformed to ROUT - R inverse 


using subroutine RINCON with dimension N + IS . If IS = 1 the subroutine 

sets RIN((N+1) (N+2))/2) = -1, so that the N+lst column of ROUT will be 

-1 1/2 

the X estimate followed by -1. R = UD so that the diagonals 
are square root scaled U columns. This information is used to con- 
struct the U-D array which is written in ROUT. 

If N<0 the input is assumed to be a U-D array. This array is 

1/2 

converted to ROUT = UD and then using RINCON, R is computed and stored 
in ROUT. If IS = 1 the U-D matrix is assumed augmented by X (estimate), 
and on output the right side term of the SRIF array is obtained. When 
IS = 1, the initial value of RIN((N+1) (N+2) /2) is restored before exiting 
the subroutine. 
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19. SFU (Sparse F * unit upper triangular U) 

Purpose 

To efficiently form the product F*U so that only the nonzero 
elements of F are employed and so that the structure of the U 
matrix is utilized (upper triangular with implicit unit diag- 
onal elements) . When F is sparse there are significant savings 
in storage and computaton. Note that since we deal only with 
the nonzero elements of F we are saved the time associated with 
computing unnecessary F matrix element addresses. 


CALL SFU (FEL , IROW , JCOL , NF , U , N , FU , MAXFU , IFU , JD I AG) 


Argument Definitions 


FEL CNF) 

IROW CNF) 
JCOL (NF) 

NF 

U(N* (N+l) / 2) 
N 

FU (IFU,N) 

MAXFU 

IFU 


JDIAG (N) 


Values of the non-zero elements of the F matrix 
Row indices of the F elements 
Column indices of the F elements 
F (IROW (K) , JCOL (IQ) - FEL (K) 

The number of non-zero elements of the F matrix 

Upper triangular, vector stored matrix with 
implicity defined unit diagonal elements. Note 
that UCJJi terms are not, in fact, unity. 

Dimension of the U matrix 

The output result 

Row dimension of the FU matrix 

Number of rows in FU. IFU. LE. MAXFU, and IFU.GE. 

Max (IROW(K), K=1,...,NF); i.e, FU must have at 
least as many rows as does F. Additional rows of 
FU could correspond to .zero rows of F. 

Diagonal element indices of a vector stored upper 
triangular matrix, i.e. JDIAG (K)=K* (JC+1)/2=JDIAG(K-1)+K. 
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ORIGINAL PAGE B 

OF POOR QUALITY 


Example: 

F(3,12) with: F(l,l) = .9, FC2,2) = .8, F(3,3) = 1.1, 
F(l,7) = 1.7, F(2,8) =-2.8 and F(3,U) = 3.11. 

In this case F has NF = 6 (nonzero elements) ; and one may 


take 


IROW(l) = 1 JCOL(l) 

IROW(2) = 2 JCOL(2) 

IR0W(3) = 3 JC0L(3) 

IR0W(4) = 1 JC0L(4) 

IR0W(5) = 2 JCOL(5) 

IR0W(6) = 3 JC0L(6) 


1 

FEL(l) 

= .9 

2 

FEL(2) 

= .8 

3 

FEL(3) 

= 1.1 

7 

FEL(4) 

= 1.7 

8 

FEL(5) 

=-2.8 

11 

FEL(6) 

= 3.11 


Remarks and'Restrictions 


Comments regarding increased efficiency are included in the code. 


Functional Description 


We write 
F 


E 


T 

F^. e. et 
ij i 1 


where e. is the i-th unit vector. Then 

X 

FU = Y) F . . e . (e?U) 

The code is based on this equation. 
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20. TDHHT (Two dimensional Householder triangularization) 

Purpose 

To transform a two dimensional rectangular matrix to a 
triangular, or partially triangular form by Householder orthogonal 
matrix pre-multiplication. This subroutine can be used to compress 
overdetermined linear systems to triangular (double subscripted 
form) in much the same way as does the subroutine THH (which outputs 
a vector subscripted triangular result) . For recursive applications 
THH is computationally more efficient and requires less storage. 

The chief application, that we have in mind, for this subroutine 
is to the matrix triangularization of "mapped" square root 
information matrices of the form S(m,n) with m less than n. 


CALL TDHHT (S ,MAXS , IRS , JCS , JSTART , JSTOP , V) 


Argument Definitions 


S (IRS, JCS) 


MAKS 


Input (possibly partially) triangular 
matrix. The output (possibly partially) 
triangular result overwrites the input . 

Row dimension of S matrix 


IRS 


Number of rows in S (IRS.LE.MAXS) , and 
IRS.GE.2. 


JCS 


Number of columns in S 


JSTART 


JSTOP 


V (IRS) 


Index of first column to be triangularized. 
If JSTART.LT. 1 then it is assumed that the 
triangularization starts at column 1. 

Index of last column to be triangularized. 
When JSTOP is not between max(l, JSTART) 
and JCS then the triangularization is 
carried out as far as possible (i.e. to IRS 
if S has less rows than columns, or to JCS 
if it has more rows than columns) . 

Work vector 
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Remarks and Restrictions 

The indices JSTART and JSTOP are input for efficiency purposes. 
When it is known that the input matrix is partially triangular one 
can by-pass the corresponding (initial) Householder reduction steps. 
Further, for certain applications it is not necessary to totally 
triangularize the input array. For example if S(m,n) and m is 
less than n, the system is in triangular form after only m elementary 
Householder reduction steps, i.e 



n 


The code is set up so that it defaults to the largest possible 


upper triangularization. 



"JSTART 

Input S 


The dotted portion of the matrix and the block of zeros are not 
employed at all in the computations. The input matrix is trans- 
formed to (possibly partially) triangular form by premultiplication 
by a sequence of elementary Householder orthogonal transformations . 


JCS 



JSTOP 
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The method is described fully in the bjooks by Lawson and Hanson 
Solving Least Squares Problems, and in Bierman - Factorization 
Methods for Discrete Sequential Estimation. 
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21 • THH 


(Triangular Householder Orthogonalization) 


Purpose 


To compute [r:z] 


such that 


T 


R 

A 


z 

z 



T - orthogonal 


. This is the key algorithm used in the square root infor ma tion batch 
sequential filter. 


CALL THH (R,N, A, IA,M, RSOS ,NSTRT) 


Argument Definitions 

R(N*(N+3)/2) Input upper triangular vector stored 

square root information matrix. If 
estimates are involved RSOS.GE.O and R 
is augmented with the right hand side 
(stored in the last N locations of R) . 

If RSOS. LT.O only the first N*(N+l)/2 
locations of R are used. The result 
of the subroutine overwrites the input R 

N Number of parameters 


A(M,N+1) Input measurement matrix. The N+lst 

column is only used if RSOS.GE.O, in 
which case it represents the right side 
of the equation v + AX = z. A is 
destroyed by the algorithm, but it is 
not explicitly set to zero. 

IA Row dimension of A 


M 

RSOS 


The number of rows of A that are to be 
combined with R (M.LE.IA) 


Accumulated residual root sum of squares 
corresponding to the data processed 
prior to this time. On exit RSOS repre- 
sents the updated root sum of squares 


of the residuals 


£ II 'i- Vest H 1 


1/2 


summed over the old and new data, 
also includes the a priori term 




It 
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2 

11 R X - z 1| . Because RSOS cannot 
1 o est o n 

be used if data, z, is not included 
we use RSOS.LT.O to indicate when data 
is not included. 

NSTART First column of the input A matrix 

that has a nonzero entry. In certain 
problems, especially those involving 
the inclusion of a priori statistics, 
it is known that the first NSTRT-1 
columns of A all have zero entries. 
This knowledge can be used to reduce 
computation. If nothing is known 
about A, then NSTRT.LE.l gives a 
default value of 1, i.e. it is assumed 
that A may have nonzero entries in the 
very first column. 


Remarks and Restrictions 

It is trivial to arrange the code so that R output need not over- 
write the input R. This was not done because, in the author's opinion, 
there are too few times when one desires to have ROUT ^ RIN . 


Functional Description 

Assume for simplicity that NSTRT=1. Then at step j, j = l,...,N 
(or N+l if data is present) the algorithm implicitly determines an 
elementary Householder orthogonal transformation which updates row j 
of R and all the columns of A to the right of the jth. At the 
completion of this step column j of A is in theory zero, but it is 
not explicitly set to zero. The orthogonalization process is discussed 
at length in the books by Lawson and Hanson - Solving Least Squares 
Problems and Bierman — Factorization Methods for Discrete Sequential 
Estimation. 
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22. TTHH 


(Two triangular matrix Householder reduction) 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Purpose 

To combine two vector stored upper triangular matrices, R and RA 
by applying Householder orthogonal transformations. The result over- 
writes R. 




j CALL TTHH (R,RA,N) 


Argument Definitions 


R(N*(N+l)/2) 


Input vector stored upper triangular 
matrix, which also houses the result 


RA(N*(N+l)/2) 


N 


Remarks and Restrictions 


Second input vector stored upper 
triangular matrix. This matrix is 
destroyed by the computation. 

Matrix dimension 

N less than zero is used to indicate 
that R and RA have right sides 
(|N|+1 columns) and have dimension 
lN|*(|N[+3)/2). 


RA is theoretically zero on output, but is not set to zero. 
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23 . TWOMAT (Triangular matrix print) 

Purpose 

To display a vector upper triangular matrix in a two 
dimensional triangular format. Precision output corresponds to a 
7 column 8 digit, double precision format. Compact output corres- 
ponds to a 12 column, 5 digit single precision format. 


CALL TWOMAT (A , N , LEN , CAR , TEXT , NCHAR , NAME S ) 


Argument Definitions 


A(N*N+l)/2) 


Vector stored upper triangular matrix (DP) 


N 


Dimension of A 


LEN Column format (7 or 12 columns) . When LEN 

is different from 7 or 12 the print defaults 
to 12 columns. 


CAR(N) 


TEXT (NCHAR) 
NCHAR 


NAMES 


Parameter names (alphanumeric) associated 
with A. When NAMES is false, CAR is not 
used. 

An array of field data characters to be 
printed as a title preceding the matrix 

Number of characters (including spaces) that 
are to be printed in text( ) 

ABS (NCHAR) .LE. 114. If NCHAR is negative there is 
no page eject before printing. NCHAR positive 
results in a page eject so that the print 
starts on a fresh page. 

A logical flag. If true then the names of 
the parameters are used as labels for the 
rows and columns. If false the output labels 
default to numerical values. 


Remarks and Restrictions 


Using NCHAR nonnegative, and starting the print at the top of a 
new page makes it easier to locate the printed result and is 
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especially recommended when dealing with large dimensioned arrays. 
Page economy can, however, be achieved using the NCHAR negative 
option. In this case the print begins on the next line. The 
alphanumerics in this routine make it machine dependent; it is 
arranged for implementation on a UNIVAC 1108. 
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24. TZERO (Triangular matrix zero) 

Purpose 

To zero out rows IS(Istart) to IF(Ifinal) of the. vector- stored 
upper triangular matrix R. 

CALL TZEROXR,N,IS,IF) 


Argument Definition 
R(N*(N+1) 12 ) 

N 

IS 

IF 

Functional Description 



R(input) 


Input vector stored upper triangular 
matrix 

Row dimension of vector stored matrix 
First row of R that is to be set to zero 
Last row of R that is to be set to zero 


IS 
IF 

R( output) 
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25. UDCOL (U-D covariance factor colored noise time update) 
Purpose 

To time update the U-D covariance factors so as to include 
the effects of colored noise variables . 


CALL UDCOL (U,N,KS , NCOLOR, V, EM, Q) 


Argument Definitions 

Input vector stored U-D covariance factors. 
The updated result resides here on output. 

Filter matrix dimension. If the last column 
of U houses the filter estimates, then 
N = number filter variables + 1. 

Location of the first colored noise variable 
(KS . GE . 1 .AND . KS -LE . N) 

The number of colored noise variables 
contiguous to the first, including the 
first. (NC0L0R.GE.1) 

Work vector ( (KS-1+NC0L0R) .LE.N) 

Input vector of colored noise mapping terms 
(unaltered by program) 

Input vector of process noise variances 
(unaltered by program) 

< 

Remarks and Restrictions 

When estimates are involved they are appended as an additional 
column to the U-D matrix. When the subroutine is applied to the 
augmented matrix the estimates are correctly updated. When the 
colored noise terms are not contiguously located one can fill in 
the gaps with unit EM terms and corresponding zero Q elements. 

It is preferable, however, to apply the subroutine repeatedly to 
the individual contiguous groups. 


U (N* (N+D/2) 

N 

KS 

NCOLOR 

> V(KS-1+NCOLOR) 
EM (NCOLOR) 

Q (NCOLOR) 
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Functional Description 

The model equation corresponding to the time update of this 
subroutine is 


X 


I 

0 

0 


X 


0 

p 

= 

0 

M 

0 


P 


i 

y 

A 

0 

0 

I 


y 


0 


L_ Jj+l 1- _J 1_ j L_ _1 

where M is diagonal, with NP terms, and w^ €N(0,Q) where Q is 
diagonal with NP terms. The output U-D array associated with this 
time update equation satisfies 

UDU T (output) = $ UDuV + BQB T 
where $ and B are as above. The algorithm for obtaining U-D 
(output) is the Bierman-Thornton one-component-at-a-time update 
described in Bierman - Factorization Methods for Discrete 
Sequential Estimation", Academic Press (1977), pp -147-148. 
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26. UDMEAS (u-D measurement update) 

Purpose 

Kalman filter measurement updating using Bierman's U-D measure- 
ment update algorithm, cf 1975 CONF. DEC. CONTROL paper. A scalar 
T 

measurement z = A x + v is processed, the covariance U-D factors 
and estimate (when included) are updated, and the Kalman gain and 
innovations variance are computed. 


CALL UDMEAS(U,N,R,A,F,G,ALPHA) 


Argument Definitions 


INPUTS 

U(N*(N+l)/2) Upper triangular vector stored input matrix. 

D elements are stored on the diagonal. The 
U vector corresponds to an a priori covariance. 

If state estimates are involved the last column 
of U contains X. In this case Dim U = (N+l)*(N+2) /2 
and on output (U(N+l)*(N+2)/2= z-A**T*X(a priori est) . 

N Dimension of state vector, N.GE.2 

R Measurement variance 


A(N) 


Vector of Measurement coefficients; if data 
then A(N+1) = z 


F(N) 


Input work vector. To economize on storage P 
can overwrite A 


ALPHA 


If ALPHA. LT. zero no estimates are computed 
(and X and z need not be included) . 


OUTPUTS 

U Updated vector stored U-D factors. When 

ALPHA (input) is nonnegative the (N+l)st 
column contains the updated estimate and 
the predicted residual. 

ALPHA Innovations variance of the measurement 

residual . 

F Contains U**T*A (input) and when ALPHA (input) 

is nonnegative F(N+1) =(z-A**T*X(a priori est))/ALPHA. 
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G(N) Vector of unweighted Kalman gains, 

K = G/ALPHA 

Remarks and Restrictions 

One can use this algorithm with R negative to delete a 
previously processed data point. One should, however, note that 
data deletion is numerically unstable and sometimes introduces 
numerical errors. 

The algorithms holds for R = 0 (a perfect measurement) and 
the code has been arranged to include this case. Such situations 
arise when there are linear constraints and in the generation of 
certain error "budgets". 

Functional Description 

The algorithm updates the columns of the U-D matrix, from 
left to right, using Bierman's algorithm, see Bierman's 
"Factorization Methods for Discrete Sequential Estimation," 
Academic Press (1977) pp 76-81 and 100-101. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


27 . UD2C0V (U-D factor to covariance) 

Purpose 

To obtain a covariance from its TJ-D factorization. Both matrices 

are vector stored and the output covariance can overwrite the input 

T 

U-D array. U-D and P are related via P = UDU . 

CALL UD2C0V (UIN, POUT ,N) 

Argument Definitions 

UIN (N* (N+l) /2) Input vector stored U-D factors, with D 

entries stored on the diagonal. 

POUT (N* (N+l) /2) Output vector stored covariance matrix 

(POUT ** UIN is permitted) . 

N Dimension of the matrices involved (N.GE.2) 
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28. UD2SIG 


(U-D factors to sigmas) 


Purpose 

To compute variances from the U-D factors of a matrix. 


CALL UD2SIG(U,N,SIG,TEXT,NCT) 


Argument Definitions 

U(N*(N+l)/2) Input vector stored array containing 

the U-D factors. The D (diagonal) 
elements are stored on the diagonal 
of U. 

N Dimension of the U matrix (N.GE.2) 

SIG(N) Output vector of standard deviations 

TEXT ( ) Output label of field data characters, 

which precedes the printed vector of 
standard deviations . 

UCT Number of characters of text, 

0.LE.NCT.LE.126. If NCT = 0, no 
sigmas are printed, i.e. nothing is 
printed. 


Remarks and Restrictions 


The user is cautioned that the text related portion of this subroutine 
may not be compatible with other computers. The changes that may be 
involved are, however, very modest. 


Functional Description 

If U and D are represented as doubly subscripted matrices then 

/ N 

SIG(J) = (D(J,J) + D(K,K) [U(J,K)] 2 

' K=J+1 

If NCT.GT.O a title is printed, followed by the sigmas. 
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29. UTINV (Upper triangular matrix inverse) 

Purpose 

To invert an upper triangular vector stored matrix and store 
the result in vector form. The algorithm is so arranged that the 
result can overwrite the input. 

CALL UTINV (RIN, N, ROUT 

Argument Definitions 

RIN(N*(N+l)/2) Input vector stored upper triangular 

matrix 

N Matrix dimension 

ROUT (N* (N+l) /2) Output vector stored upper triangular 

matrix inverse (ROUT = RIN is permitted) 

Remarks and Restrictions 

111 conditioning is not tested, but for nonsingular systems the 
result is as accurate as is the full rank Euclidean scaled 
singular value decompostiion inverse. Singularity occurs if a 
diagonal is zero . The subroutine terminates when it reaches a 
zero diagonal. The columns to the left of the zero diagonal are, 
however, inverted and the result stored in ROUT. 

This routine can also' be used to produce the solution to RX = 2. 
Place Z in column N+l (viz. RIN (N* (N+l) /2+1) = Z(l), etc.), define 
RIN((N+1) (N+2)/2) = -1 and call the subroutine using N+l instead 
of N. On return the first N entries of column N+l contain the 
solution (e.g. R0UT(N*(N+1) /2+1) = X(l), etc.). When only the 
estimate is needed, then it is more efficient to use the code 
described in section to II. 8 to obtain X, directly. 
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Because matrix inversion is numerically sensitive we recommend 
using this subroutine only in double precision. 

Functional Description 

The matrix inversion is accomplished using the standard back 
substitution method for inverting triangular matrices, cf . the book 
references by Lawson and Hanson, [1] or Bierman [3]. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


30. UTIROW (Upper triangular inverse, inverting only the upper rows) 
Purpose 

To compute the inverse of a vector stored upper triangular 
matrix, when the lower right comer triangular inverse is given. 


CALL UTIROW (RIN, N, ROUT, NRY) 


Argument Definitions 

RIN (N* (N+l) /2 ) Input vector stored upper triangular 

matrix. Only the first N - NRY rows 
are altered by the algorithm. 

N Matrix, dimension. 


ROUT (N* (N+l) /2) Output vector stored upper triangular 

matrix inverse. On input the lower 
NRY dimensional right comer contains 
the given (known) inverse. This lower 
right corner matrix is left unchanged. 
(ROUT = RIN is permitted.) 

NRY Number of rows, starting at the bottom, 

that are assumed already inverted. 

Remarks and Restrictions 


The purpose of this subroutine is to complete the computation 


of an upper triangular matrix inverse, given that the lower right 


comer has already been inverted. Part of the input, the rows to 
be inverted, are inserted via the matrix RIN. The portion of the 


matrix that has already been inverted is entered yia the matrix ROUT. 
It may seem odd that part of the input matrix is put into RIN and 

part into ROUT, The reasoning behind this decision is that RIN 

represents the input matrix to be inverted (it just happens that 


we do not make use of the lower right triangular entries); ROUT 


represents the inversion result, and therefore that portion of the 
inversion that is given should be entered in this array. 
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Ill conditioning is not tested, but for nonsingular systems the 

result is accurate. Singularity halts the algorithm i-f any of the 

first N-NRY diagonal elements is zero. If the first zero encountered 

moving up the diagonal (starting at N-NRY) is at diagonal j then the 

rows below this element will be correctly represented in ROUT. 

To generate estimates do the following: put N+l into the matrix 

dimension argument; in the first N-NRY rows of the last column of 

RIN put the right hand side elements of the equation R^x + R y = 

(i.e., R , R , and z make up the first N-NRY rows of RIN); in the 
x xy x 

next NRY entries of ROUT, beginning in the (N-NRY+l)st element, put 

y (i.e., R ^ and y make up rows N-NRY+1, , , . ,N of ROUT); and 
6S u y 6St 

ROUT ( (N+l) (N+2)/2) = -1. On output, the last column of ROUT will 

contain x _ , y and -1 . 
est’ J est 

When NRY = 0 this algorithm is equivalent to subroutine UTINV. 
Functional Description 

The matrix inversion is accomplished using the standard back 
substitution method. The computations are arranged, row-wise, starting 
at the bottom (from row N-NRY, since it is assumed that the last NRY 
rows have already been inverted) . 
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31. WGS (Weighted Gram— Schmidt matrix triangularization) 

Purpose 

To compute a vector stored U-D array from an input rectangular 

T T 

matrix W, and a diagonal matrix D so that W D W = UDU . 
56 w w 

GALL WGS (W, IMAXW, IW,JW,DW,U,V) 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Argument Definitions 
W(IW,JW) 

IMAXW 

IW 

JW 

DW(JW) 


U(IW*(IWfl)/2) 

V(JW) 

Remarks and Restrictions 


Input rectangular matrix, destroyed by 
the computations 

Row dimension of input W matrix, 

IMAXW. GE.IW 

Number of rows of W matrix, dimension of U 

Number of columns of U matrix 

Diagonal input matrix; the entries are 
assumed to be nonnegative. This vector 
is unaltered by the computations 

Vector stored output U-D array 

Work vector in the computation 


The algorithm is not numerically stable when negative DW weights 

are used; negative weights ar§, however, allowed. If JW is less than 
IW (more rows than columns), the output U-D array is singular; with 
IW-JW zero diagonal entries in the output U array. 

Functional Description 

A D^-orthogonal set of row vectors, <j>^, ^xw* are con ~ 

T 

structed from the input rows of the W matrix, i.e., W = U <j>, , <j>D <j> = D. 

w 

The construction is accomplished using the modified Gram-Schmidt 
orthogonal construction (see refs. [1] or [3]). This algorithm is 

reputed to have excellent numerical properties. Note that the <f> 

vectors are not of interest in this routine, and they are overwritten; 

The V vector used in the program houses vector IW-j+1 of 4> at step j of 

algorithm. The fact that the computed <)> vectors may not be D orthogonal 
is of no import in regard to the U and D computed results. 
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V. FORTRAN Subroutine Listings 

The subroutines use only FORTRAN IV, and are therefore essentially 
portable. The one notable exception is subroutine TWOMAT, which prints 
triangular, vector stored matrices. It employs FORTRAN V FORMAT state- 
ments and six character UNIVAC alphanumeric wordlength, and thus is UNIVAC 
dependent. Subroutine TJD2SIG also involves text, and it too is therefore 
to some extent machine dependent. Comment statements appear occasionally 
to the right of the FORTRAN code, and are preceded by a symbol. The 
subroutine user can, if necessary, transfer or remove such program 
commentary. 

All of the subroutines employ "implicit double precision” statements . 
They are, however, constructed so as to operate in single precision, and 
the user has only to omit or comment out the implicit statements. If the 
subroutines are to be used in double precision on a machine that does not 
have the implicit FORTRAN option one should explicitly declare all of the 
non-integer variable names appearing in the programs as double precision 
variables . 

If these subroutines are to be used in production code and computa- 
tional efficiency is of major concern one should replace the somewhat 
lengthy subroutine argument lists by introducing COMMON, and including 
those terms in the COMMON that are redundantly computed with each sub- 
routine call. 
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ooooooooooooooooonooo 


SUBROUTINE A2A1 ( A . 1 A . IR.LA.NAMA. Ai . IA1 .LAI .NAMA1 ) 

SUBROUTINE TO REARRANGE THE COLUMNS OF A(IR.LA). IN NAMA ORDER 
AND PUT THE RESULT IN Al(IR.LAl) IN NAMA1 ORDER. ZERO COLUMNS 
ARE INSERTED IN Al CORRESPONDING TO THE NEWLY DEFINED NAMpS, 


A(lRfLA) 

IA 

IR 

LA 

NAM/\ (LA) 
AKlR.LAl) 

IA1 

LAI 

NAMfll(LAl) 


input rectangular MATRIX 

ROW DIMENSION of a. IR.LE.IA 

MO. OF ROWS OF A THAT ARE TO BE REARRANGED 

NO. COLUMNS IN A * ALSO THE 

NO. OF PARAMETER NAMES ASSOCIATED WITH A 

PARAMFTER NAMES ASSOCIATED WITH A 

output rectangular matrix 

A AND Al CANNOT SHARE COMMON STORAGE 

ROW DIMENSION OF Al. iR.LE.IAl 

NO. COLUMNS IN Al. ALSO THE 

NO. OF PARAMETER NAMES ASSOCIATED WITH Al 

INPUT LIST OF PARAMETER NAMES TO BE ASSOCIATED 

WITH THE OUTPUT MATRIX Al 


COGNIZANT PERSONS! G, J.BIERMAN/M« W.NEAO (JPL» SEPT. 1976) 

DIMENSION A( IA. 1 ) » MAMA ( 1 ) » Al < IA1 . 1 > »NAMA1 ( 1 ) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

ZERO=0. 

DO 100 J=i.LAl 
DO 6g 1=1. LA 

IF (NAMA(I) .EQ.nAMAKJ)) GO TO 80 
60 CONTINUE 

DO 7q K=1 ► IR 
70 Al(K»J)=ZERO 

GO TO 100 
80 DO 9f) K=1 » IR 
90 A1(K.J)=A(K.I) 

100 CONTINUE 

RETURN 
END 


0 ZERO COL. CORRES. TO NFW NAME 
9 COPY COL. ASSOC. WITH OLD NAME 


A2A1001C 

A2A10020 

A2A10030 

A2A10040 

A2A10050 

A2A10060 

A2A10070 

A2A100B0 

A2A10090 

A2A10100 

A2A10110 

A2A10120 

A2A10130 

A2A10140 

A2A10150 

A2A10160 

A2A10170 

A2A10180 

A2A10190 

A2A10200 

A2A10210 

A2A10220 

A2A10230 

A2A10240 

A2A10250 

A2A10260 

A2A10270 

A2A10280 

A2A10290 

A2A10300 

A2A10310 

A2A10320 

A2A10330 

A2A10340 

A2A10350 

A2A10360 

A2A10370 

A2A10380 

A2A10390 
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SUBROUTINE COMBO (R »Ll »NAM1 ,L2 ,NAM2, A# I A *LA »NAMA ) 


COMBOnOO 


C 


c 


50 


60 


70 

80 

100 


150 


160 


TO REARRANGE A VECTOR STORED TRIANGULAR MATRIX AND STORE 

the result in matrix a. the difference between this sub- 
routine AND R2A IS THAT THERE THE NAMELTST ^0R A IS INPUT, 
HERE IT IS DETERMINED BY COMBINING THE LIST FOP R WITH 
A LIST OF DESIRED NAMES. 


R<H*<Ll+l)/2) 

LI 

NAMl(Ll) 

L2 

NAM2(L2> 


A (Ll » LA ) 

IA 

LA 


NAMa(LA) 


INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
NO. OF PARAMETERS IN R t AND IN NAMl ) 

NAMES ASSOCIATED WITH R 
NO. OF PARAMETERS IN NAM2 

parameter names that ARE to BE combined with 

(NAMl LIST). THESE NAMES MAY OR MAY NOT RE IN 

n ami . 

output ARRAY containims the rearranged 
r matrix* li.le.ia. 

ROW DIMENSION of a 

NO. OF PARAMETER NAMES IN NAMA* AND THE 
COLUMN DIMENSION Of A. LA=L1+L2-N0, nam^s 
COMMON TO NAMl AnD NAM2. LA IS COMPUTED AND 
OUTPUT. 

PARAMETER NAMES ASSOCIATED with THE OUTPUT A 
MATRIX. CONSISTS OF NAMES IN NAMl WHICH ARE 
NOT IN NAM2 FOLLOWED BY NAM2. 


COGNIZANT PERSONS? G.J.BIERMAN/M.W.NEAD ( JPL * SEPT. 1976) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 

DIMENSION R(ll> A(IA*1>» NAM1(1)» nAM 2(1). NAMA(l) 


ZERO=0.0 

K=1 

DO 100 I=1*L1 
DO 5q J=1*L2 

IF (NAMl(I) ,EQ.NAM2(J) ) GO TO 1Q0 
CONTINUE 
NAMA(K)=NAM1(I) 

JJ=I*(I-l)/2 
DO 6o L=l* I 

A(L*K)=R(JJ+U 
IF (I.E0.L1) GO To 80 
IP1 = 1+1 
DO 70 L=IP1*L1 
A(L*K> = ZERO 
K=K+1 
CONTINUE 

NAMES UNIQUE TO NAMl ARE NOW In NAMA 
DO 200 J=1»L2 
DO 150 I±1»L1 

IF (NAM2(J) .EQ.NAMl(D) GO TO 170 
CONTINUE 
NAMA(K)=NAM2(J) 

DO 160 L=1*L1 
A(L'K)=ZERO 


COMBOOIO 
COMBO020 
COMB0030 
COMBOP40 
COMBO05O 
COMB0060 
COMB0070 
COMB0080 
COMBO090 
COMBOIOO 
R COMFOllO 
COMB0120 
COMBO) 30 
COMB0140 
COMBO150 
COMB0160 
COMB0170 
COM0O18O 
COMBO190 
COMB0200 
COMBO210 
COMB0220 
COMB0230 
C0MB0240 
COMBO250 
COMBO260 
C0MR0270 
COMBO280 
COMB0290 
COMB0300 
C0MB0310 
COMB0320 
COMB0330 
COMB0340 
COMB0350 
COMBO360 
C0MB0370 
COMBO380 
COMBO390 
COMBO400 
COMB0410 
COMBO420 
COMB0430 
COMBO440 
COMBO450 
C0MBO460 
COMB0470 
COMB0480 
COMBO490 
COMBO500 
COMBO510 
COMBO520 
COMB0530 
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C NAMES UNIQUE TO NAM? ARE NOW IN NAMA 

GO TO 190 

170 NAMA(K)=NAM2(J> 

C LOCATE DIAGONAL OF PRECEDING COLUMN 

JJ=I*<I~l)/2 > 

DO 180 L=lfl 
180 A(L»K)=R(JJ+L) 

IF (I.EQ»LD GO TO J.90 
IP1=I+1 

DO 185 L=IP1»L1 
185 A(L*K)=ZERO 
190 K=K+1 

200 CONTINUE 
LA=K-1 

C ' NAMES MUTUAL TO NAMl AND NAM2 ARE NOW IN NAMA 

RETURN 
• END 


COMB0540 

COMB0550 

COMB0560 

COMBO570 

COMBO580 

COMB0590 

COMB0600 

COMSO610 

COMB0620 

COMBO630 

COMBO640 

COMBO650 

COMB0660 

COMBO670 

COMB0680, 

COMBO690 

COMBOTOO 


88 



ooo o n onnooooooooooo 


ORIGINAL PAGE IS 
OF POOR QUALITY 


10 


20 


SUBROUTINE C0VRH0 <C0VrN>RH0»V) 

TO COMPUTE THE CORRELATION MATRIX RH0» FROM AN INPUT COVARIANCE 
MATRIX COV. BOTH MATRICES ARE UPPER TRIANGULAR VECTOR STORED, 

The CORRELATION MATRIX result CAN OVERWRITE the INPUT COVAPIANCE 
C0V(N*(N+l>/2) INPUT VECTOR STORED POSITIVE SEMI-DFFINITE 
COVARIANCE MATRIX 

N NUMBER OF PARAMETERS r N.GE.l 

RH0(N<N+i)/2> OUTPUT VECTOR STORED CORRELATION MATRIX t 
RH0(IJ)=C0V(I<J3/(SIGMA(I)*SIGMA(J) ) 

V(N) WORK VECTOR 

COGNIZANT PERSONS: G.J.BIERMAN/M*W.NEAD t JPLfFEB.1978) 

IMPLICIT DOUBLE PRECISION <A-H»0-Z) 

DIMENSION C0V(1>* RH0(1)» Vd) 

ONE=1.DO 

Z=O.D0 

JJ=0 

DO 10 J=1»N 
JJSJJ+J 
V(J)=Z 

IF (COV(JJ).GT.Z) V(J)=0NE/ SORT (COV < JJ> ) 

**** some machines require dsqrt for double precision 

CONTINUE 
I J=0 

DO 20 J— 1 * N 
S=V(J) 

DO 20 I=1*J 
IJ=IU+1 

RHO ( I J ) = C0V ( I J) *S*V ( I ) 

RETURN 

END 


COVRHOIO 

COVRH020 

COVRH030 

COVRH040 

COVRH050 

COVRH060 

COVRH070 

COVRH080 

COVRH090 

COVRHIOO 

COVRH110 

C0VRH120 

COVRH130 

COVRH140 

COVRH1SO 

C0VRH160 

COVRH170 

COVRH180 

COVRH190 

COVRH20O 

COVRH210 

COVRH220 

C0VRH230 

COVRH240 

COVRH250 

C0VRH260 

C0VRH270 

COVRH280 

COVRH290 

COVRH300 

COVRH310 

COVRH320 

COVRH330 

C0VRH340 

COVRH350 

COVRH360 

C0VRH370 

COVRH380 

COVRH390 
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SUBROUTINE C0V2RI (U»N) 

. TO CONSTRUCT THE UPPER TRIANGULAR CHOLESKY FACTOR OF A 
POSITIVE SEMI-DEFINITE MATRIX. BOTH THE INPUT COVARIANCE 
AND THE OUTPUT CHOLESKY FACTOR (SQUARE ROOT) ARE VECTOR 
STOREO. THE OUTPUT OVERWRITES THE INPUT. 
COVARIANCE(INPUT)=U*U**T (U IS OUTPUT). 

IF THE INPUT COVARIANCE IS SINGULAR THE OUTPUT FACTOR HAS 
ZERO COLUMNS. 

U(N*(N+l>/2) CONTAINS THE INPUT VECTOR STORED COVARIANCE 

matrix (assumed positive definite) and on output 

IT CONTAINS the upper TRIANGULAR SQUARE ROOT 
FACTOR. 

, N DIMENSION OF THE MATRICES INVOLVED 

COGNIZANT PERSONS*. G.J.BIERMAN/M.W.NEAD (JPL> FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION U(l) 

ZERO=0.0 

ONE=l. 

JJ=N*(N+l)/2 
DO 5 J=N»2»-1 

IF (U(JJ) .LT»ZERO) U(JJ)=ZEPO 
U(JJ)= SORT (U(Jj) ) 

ALPHA=ZER0 

IF (U(JJ) .GT.ZERO) ALPHA=ONE/U(jJ) 


0 JJN+K=(K»J) 


KK=0 

JJN=JJ-J 0 NEXT DIAGONAL 

JM1=U-1 

DO 4 K=1»JM1 

UUJN+K)=ALPHA*U(JJN+K) 0 JJN+K=(K»J 

S=U(JJN+K) 

DO 3 I=1»K 

3- U(KK+I)=U(KK+I)-S*IJ(JJN+I) Q KK+I=(T.K) 

4 KK=KK+K 

5 UJ=JJN 

IF (U(l).LT.ZERO) UCl)=ZERO 
U(l)= SQRT(U(1) ) 

RETURN 
ENO 


COV2R010 

COV2R020 

COV2R030 

COV2R040 

COV2RO50 

COV2R060 

COV2R070 

COV2R080 

COV2R090 

COV2R100 

COV2R110 

COV2R120 

COV2R130 

COV2R140 

COV2R150 

COV2R160 

COV2R170 

C0V2R180 

COV2R190 

COV2R200 

COV2R210 

C0V2R220 

C0V2R230 

C0V2R240 

COV2R250 

COV2R260 

COV2R270 

COV2R280 

COV2R290 

CCV2R300 

COV2R310 

COV2R320 

COV2R330 

COV2R340 

C0V2R350 

COV2R360 

C0V2R370 

COV2R380 

COV2R390 

COV2R400 

COV2R410 

COV2R420 

COV2R430 

COV2R440 

COV2R450 

COV2R460 

COV2R470 
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C 

C 

C 


SUBROUTINE C0V2U0 <U#N) COV2U010 

COV2U020 

TO OBTAIN THE U-D FACTORS OF A POSITIVE SEMI-DEFINITE MATRIX. C0V2U030 
THE INPUT VECTOR STORED MATRIX IS OVERWRITTEN BY THE OUTPUT COV2U040 
U-D FACTORS WHICH ARE ALSO VECTOR STORED. COV2U050 

COV2U060 

U(N*(N+l>/2> CONTAINS INPUT VECTOR STORED COVARIANCE MATRIX. COV2U070 

ON OUTPUT IT CONTAINS THE VECTOR STORED U-D COV2UO80 

COVARIANCE FACTORS. COV2U090 

N MATRIX DIMENSION# N#GE.2 COV2U100 

COV2U110 

SINGULAR INPUT COVARIANCES RESULT IN OUTPUT MATRICES WITH ZERO COV2U120 
COLUMNS COV2U130 

COV2U140 

COV2U150 

COGNIZANT PERSONS: G.J.BIERMAM/R. A. JACOBSON ( JPL# FEB. 1977) COV2UI60 

COV2U170 

IMPLICIT DOUBLE PRECISION (A-H#0-Z) COV2U180 

COV2U190 

DIMENSION U(l) COV2U2O0 

COV2U210 


Z=O.DO 

ONE=1.DO 

NONE=l 


COV2U220 

COV2U230 

COV2U240 

COV2U250 


JJ=N* ( N + l) /2 
NP2=N+2 
DO 50 L=2*N 
J=NP2-L 

alpha=z 

IF (U(JJ).GE.Z) GO TO 10 
WRITE (6#100) J#U(JJ> 

U(JJ)=Z 

10 IF (utJJ).GT.Z) ALPHA=ONE/U( JJ) 
J J—JJ—J 
KK=0 
KJ=JJ 
JM1-J— 1 
DO 4o K=1»JM1 
KJ=KJ+1 
BETA=U(KJ) 

U(KJ)=ALPHA*U(KJ> 

I Jr JJ 
IK=KK 

DO 30 I = 1»K 
IK=IK+1 
IJ=IJ+1 

30 u(ik)=u(ik)-beta*u(ij) 

HO KKrKK+K 

50 CONTINUE 

IF (U(l).GE.Z) GO TO 60 
WRITE (6#100) NONE# U(l> 

U ( 1 ) =z 

60 RETURN 


COV2U260 

C0V2U270 

COV2U280 

COV2U290 

COV2U300 

COV2U310 

COV2U320 

COV2U330 

COV2U340 

COV2U350 

COV2U360 

COV2U370 

COV2U380 

COV2U390 

COV2U400 

COV2U410 

COV2U420 

COV2U430 

COV2U440 

COV2U450 

COV2U460 

C0V2U470 

COV2U480 

C0V2U490 

COV2U500 

COV2U510 

COV2U520 

COV2U530 

COV2U540 

COV2U550 
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COV2U560 

COV2U570 


100 FORMAT ( 1H0 1 20X t * AT STEP' »l4» ’DIAGONAL ENTRY ='»E12.4) 
END 
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10 

30 


40 


80 


SUBROUTINE C2C (C > IC * LI »NAM1 *L2 • NAM2 > 

subroutine to rearrange the ROWS And columns of matrix 

C(Li*Ll) IN NAM1 ORDER AND PUT THE RESULT IN 
C(L2»L2) IN NAM2 ORDER. ZERO COLUMNS AND ROWS ARE 
ASSOCIATED WITH OUTPUT DEEINED NAMES THAT ARE NOT CONTAINED 
IN NAM1. 

C(L1»L1) INPUT MATRIX 

IC ROW DIMENSION OF C * lC.GE.L=MAX(Ll*L2> 

LI NO. OF PARAMETER NAMES ASSOCIATED WITH THE INPUT C 

NAM1(L) PARAMETER NAMES ASSOClATEO WTTH C ON INPUT. (ONLY 

THE FIRST LI ENTRIES Appj_Y TO THE INPUT C) 

L2 NO. OF PARAMETER NAMES ASSOCIATED WITH THE OUTPUT 

NAM2(L2) PARAMETER NAMES ASSOCIATED WITH THE OUTPUT C 

COGNIZANT persons: G.J.bIERMAN/M.W.NEAD (JPL> SEPT. 1976) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 

DIMENSION C(IC»l)r NAMl(l ) t NAM211) 

ZERO=0. 

L=MAX(Ll»L2) 

IF CL.LE.L1) GO TO 5 
NM=L1+1 
DO 1 K=NM»L 
NAM1(K)= ZERO 
DO 90 J=1>L2 
DO 10 1-1 » L 

IF (NAMl(I) .EQ.NAM2CJ) ) GO TO 30 
CONTINUE 


0 ZERO REMAINING NAM1 LOCNS 


90 


GO TO 90 

IF (I.EQ.J) GO TO 90 
DO 4o K=1.L 
H=C<Kr J) 
C(K»J)=C(K»I) 

c(k*i)=h 

DO 80 K=1»L 
H=C(J»K) 
C(j»K)=C(IrK) 

C ( I »K)=H 
NMrNAMl ( I ) 
NAM1(I)=NAM1 (J) 
NAMl(J)=NM 
CONTINUE 


(3 interchange columns i and j 


0 INTERCHANGE ROWS I AND U 


ra interchange labels i and j 


100 


FIND NAM2 NAMES NOT IN NAM1 AMD SET CORPESPONDING ROWS AND 

columns TO zero 

DO 120 J=1»L2 
DO loO I=1»L 

IF (NAMl(I) .EQ.NAM2(J) ) GO TO 120 
CONTINUE 
DO HO K=lfL2 
C(J»K)=ZERO 


C2COOOOO 

C2C00010 

C2C00020 

C2C00030 

C2C00040 

C2C00050 

C2C00060 

C2C00070 

C2C00080 

C2C00090 

C2C00100 

C2C00110 

C2C00120 

CC2C00130 

C2C00140 

C2C00150 

C2C00I60 

C2C00170 

C2C00180 

C2C00190 

C2C00200 

C2C00210 

C2C00220 

C2C00230 

C2C00240 

C2C00250 

C2C00260 

C2C00270 

C2C00280 

C2C00290 

C2C00300 

C2C00310 

C2C00320 

C?C00^30 

C2C00340 

C2C00350 
C2C00360 
C2C00370 
C2C00380 
C2C00390 
C2C00400 
C2C0041 0 
C2C00420 
C2C00430 
C2C00440 
C2C00450 
C2C00460 
C2C00470 
C2C00480 
C2C00490 
C2C00500 
C2COOS10 
C2C00520 
C2C0OS30 
C2C00^40 
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110 C<K'0)=ZERO C2C00550 

120 CONTINUE C2C00560 

C C2C00S70 

RETURN C2C00580 

END C2C00590 
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SUBROUTINE HHPOST ( S » W » MROW , NPOW t NCOL » V ) 

TRIAN6ULARIZES RECTANGULAR W BY POST MULTIPLYING TT BY An 
ORTHOGONAL TRANSPORMATION T. THE RESULT IS IN S 


S(NR0W*(NR0W+l)/2) 

W(NROWfNCOL) 

mrow 

NROW 

NCOL 
V (NCOL) 


OUTPUT UPPER TRIANGULAR VECTOR STORED SQRT 
COVARIANCE MATRIX 

INPUT RECTANGULAR SORT COVARIANCE MATRIX 
(W IS DESTROYED BY COMPUTATIONS) 

ROW DIMENSION OF W 

number or Rows or w to bf trtangularized 
and THE DIMENSION OF S 
NUMBER OF COLUMNS OF W 
WORK VECTOR 


(NROW.GT.l) 
(MC0L.GE.MROW) 


COGNIZANT PERSONS: G.J.BIERMAN/M.W.nEAD (JPLf NOV. 1977) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) 

DOUBLE PRECISION SUM > BETA 
DIMENSION S(l) i W (MROW * NCOL) »V(NCOL) 

ZERO=O.DO 

one=i.do 

JC0L=NC0L 

NSYM=NR0W* ( NROW+l ) /2 
JC=NR0w+2 
DO 150 L=2»NR0W 
IR0W=JC-L 
SUM=ZERO 
DO 1qO K = 1 f JCOL 
V(K)=W(IROW^K) 

100 SUM=SUM+V(K)**2 

SUM=dSQRT(SUM) 

IF (V(JCOL) .GT.ZERO) SUMs-SUM 0 DIAGONAL ENTRY (JCOL»JCOL) 
S(NSYM)=SUM 

nsym=nsym-irow 

V ( JCOL ) =V ( JCOL ) -SUM 

IF (SUM. NE. ZERO) BETA=-ONE/<SUM*V(JCOL> > 

C T(ORTHOG. trans. )=I-BETA*V*V**T 

IR0Wm 1=IR0W-1 
JCOLMl=JCOL-l 
DO 140 I = 1 > IR0WM1 
SUM=2ER0 
DO 110 K=1>JC0L 
110 SUM=SUM+V(K)*W(I»K) 

SUM=BETA*SUM 
DO 120 K=1»JC0LM1 
120 W(I»K)=W(I»K)-SUM*V(K) 

140 S ( NSYM+I ) =W ( I » IROW ) -SUM*V ( IROW) 

150 JC0L=JC0LM1 
C 

JC=NC0L-NR0W+1 

sum=zero 


HHPOSOlO 

HHPosn20 

HHPOS030 

HHPOS040 

HHPOS050 

HHPOS060 

HHPOSOTO 

HHPOS080 

HHPOS090 

HHPOS100 

HHPOS110 

HHPOS120 

HHPOSI30 

HHPOS140 

HHP05150 

HHPOS160 

HHPOS170 

HHP0S180 

HHPOS190 

HHPOS200 

HHPOS210 

HHP0S220 

HHP0S230 

HHP0S24C! 

HHPOS250 

HHPOS260 

HHPOS270 

HHPOS280 

HHPOS290 

HHP05300 

HHPOS310 

HHP0S320 

HHPOS330 

HHPOS340 

HHP0S350 

HHPOS360 

HHPOS370 

HHPOS380 

HHPOS390 

HHPOS400 

HHPOS410 

HHPOS420 

HHPOS430 

HHPOS440 

HHPOS450 

HHPOS460 

HHPOS470 

HHPOS480 

HHP0S490 

HHP05500 

HHPOS510 

HHPOS520 

HHPOS530 

HHPOS540 

HHPOS550 
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00 160 J=1fJC HHPOS560 

160 SUM=SUM+W<I»J)**2 HHPOS570 

S ( 1 ) =DsQRT (SUM) HHPOS580 

C HHPOS590 

RETURN HHPOS600 

END HHPOS610 


96 



nno o o o o o o o r> r> o o o o o o r> o o o 


ORIGINAL PAGE IS 
OP POOR QUALITY 
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SUBROUTINE INF2R (R.N) 

TO CHOLESKY FACTOR AN INFORMATION MATRIX 
COMPUTES A LOWER TRIANGULAR VECTOR STORED CHOLESKY FACTORIZATION 

of a positive semi-definite matrix. r=r(**t)r. r upper triangular. 

BOTH MATRICES ARE VECTOR STORED AND THE RESULT OVERWRITES 
THE INPUT 

R(N*(N+l)/2) ON INPUT THIS IS A POSITIVE SEMI -DEFINITE 

UNFORMATION) MATRIX. AND ON OUTPUT IT IS THE 
TRANSPOSED LOWER TRIANGULAR CHOLESKY FACTOR. IF THE 
INPUT MATRtX IS SINGULAR THE OUTPUT MATRIX WILL 
HAVE ZERO DIAGONAL ENTRIES 

N DIMENSION OF MATRICES INVOLVED » N.GE.2 

COGNIZANT PERSONJ G. J.BIERMAN/M.W.NeAD (JPL. FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION R(l) 

Z=0.D0 

ONE=l.D0 

JJ-0 

NN=N*<N+1)/2 

NMI=N-1 

DO 10 J=i»NMl 

UJ=JJ+J 0 JJ=(J.J) 

IF (R(JJ).GE.Z) GO TO 5 
WRITE (6.20) J.R(JJ) 

R(JJ)=Z 

5 R(JJ)= SQRTCR(JJ)) 

**** SOME MACHINES REQUIRE DSQRT FOR DOUBLE PRECISION 
ALPHA=Z 

IF (R(JJ).GT.Z) alpha=one/pcjj) 

JK=NN + ^ 0 jKr(j.K) 

JP1=J+1 

JIS=JK Q JI5={J»I) START 

NPUP1=N+JP1 

DO lo L=JP1»N 

k=npjpi-l 

JK=JK-K 

R(JK)=ALPHA*R(JK) 

beta=r(jk) 

ki=nn+k 

JI=JIS 
NPk=N+K 
DO 10 M=K»N 
I=NPK-M 
Kl=KI-I 


INF2R010 

INF2R020 

INF2R030 

INF2R040 

INF2R050 

INF2R060 

INF2R070 

INF2R080 

INF2R090 

INF2R1C0 

INF2RU0 

INF2R120 

INF2R130 

INF2RJ40 

INF2R150 

INF2R160 

INF2R170 

INF2R180 

IMF2R190 

INF2R200 

IMF2R210 

INF2R220 

INF2R230 

INF2R240 

INF2R250 

INF2R260 

INF2R270 

INF2R280 

INF2R290 

INF2R300 

INF2R310 

INF2R320 

INF2R330 

IMF2R340 

IMF2R350 

IMF2R360 

INF2R370 

IMF2R380 

INF2R390 

INF2R400 

IMF2R410 

INF2R420 

IMF2R430 

INF2R440 

IMF2R450 

INF2R460 

IMF2R470 

IMF2R480 

INF2R490 

IMF2RSOO 

INF2R510 

INF2R520 

IMF2R530 

INF2RS40 

INF2R550 
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10 . R-CKI-)-=R-(K-I-)-R{JlM'BETA - INF2R560 

c INF2R570 

IF (R(nN).GE.Z) GO TO 15 IMF2R580 

WRITE ( 6»20 ) NfR(NN) IMF2R590 

R<NN)=Z INF2R600 

15 R(NN)= SORT <R (NN) ) IMF2R610 

RETURN INF2R620 

C INF2R630 

20 FORMAT (1HO»20Xf’ AT STEP »» 14 »* DIAGONAL ENTRY =*»E12.4, INF2R640 

1 *» IT IS RESET TO ZERO') INF2R650 

END ' IMF2R660 
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ORIGINAL PA.GEB 
OF POOR QUALITY 


c 


c 


c 


SUBROUTINE PERMUT (A»IArIR»Ll»NAMl rL2fNAM2) 

SUBROUTINE TO REARRANGE PARAMETERS OF A(IR,L1>. MAM1 ORDER 
TO A<IR»L2)> NAM2 ORDER. ZERO COLUMNS ARE INSERTED 
CORRESPONDING TO THE NEWLY DEFINED NAMES. 

AdRrL) INPUT RECTANGULAR MATRIX. L=MAX<Ll»L2) 

I A ROW DIMENSION OF A* IA.GE.IR 

IR NUMBER OF ROWS OF A THAT ARE TO BE REARRANGED 

LI NUMBER OF PARAMETER NAMES ASSOCIATED WITH THE INPUT 

A MATRIX 

NAM1(L) PARAMETER NAMES ASSOCIATED WITH A ON INPUT 

(ONLY THE FIRST LI ENTRIES APPLY TO THE INPUT A) 

NAM1 IS DESTROYED RY PfRMUT 

L2 NUMBER OF PARAMETER NAMES ASSOCIATED WITH THE OUTPUT 

A MATRIX 

NAM2 PARAMETER NAMES ASSOCIATED WITH THE OUTPUT A 

COGNIZANT PERSONS? G, J.RlERMAM/M.W.MEAD <JPL» SEPT.- 1976) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A ( I A » 1 ) t nAMKD. NAM2U) 

ZERO=0, 

L=MAX(Ll»L2) 

IF (L.LE.Ll) GO TO 50 
NM=L1+1 
- DO 40 K r NM»L 
40 NAM1<K)=0 

50 DO 100 J=1»L2 
DO 6o 1=1 >L 

IF (NAMl(I) ,E0.NAM2(J) ) GO TO 65 
60 CONTINUE 
GO TO 100 
65 CONTINUE 

IF ( I »EQ. J) GO TO 100 

DO 7q K=1»IR R INTERCHANGE COLS I AND J 


D ZERO REMAINING NAM1 LOGS 


PFRMU010 

PERMU020 

PERMU030 

PERMU040 

PFRMU050 

PERMU060 

PERMU070 

PERMU080 

PF.RMU090 

PFRMU100 

PERMU110 

PERMU120 

PERMU130 

PERMU140 

PERMU150 

PERMU160 

‘PFRMU170 

PERMU180 

PERMU190 

PERMU200 

PFRMU210 

PERMU220 

PERMU230 

PERMU240 

PERMU250 

PERMU260 

PFRMU270 

PERMU280 

PFRMU290 

PrRMUSOO 

PERMU310 

PERMU320 

PERMU330 

PERMU340 

PERMU350 

PERMU360 

PFRMU370 


W=A(K»J) 

A(K»J>=A(K. I) 

70 A(K»I)=W 

NMrNAMl(I) fl INTERCHANGE I AND J COL. LABELS 

NAm1(I)=NAM1(J) 

- NAMl(J)=NM 
100 CONTINUE 

REPEAT TO FILL NFW COLS 

DO 2o0 J=1.L2 
DO 160 1=1 .L 

IF (NAMl(I) .E0.NAM2(J) ) GO TO 200 
160 CONTINUE 

DO 170 K=1.IR 
170 A(K»J)=ZERO 

200 CONTINUE 

RETURN 

END 


PERMU380 

PERMU390 

PERMU400 

PERMU410 

‘ PERMU420 
PERMU430 
PFRMU440 
PERMU450 
PFRMU460 
PERMU470 
PERMU480 
PERMUU90 
PERMU500 
PEPMU510 
PERMU520 
PERMU530 
PERMU540 
PERMU550 
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C 

C 


C 


■SUBROUTINE PHIUfPHI t MAXPHI t iRPHI t ICPHJ #U»N#PHIU»MPHIU> 


PHIU001O 

PHIU0020 


this subroutine computes w=phi*u whe r e phi is a rectangular matrixphiuoo30 
'WITH IMPLICITLY DEFINED COLUMNS OF TRAILING ZEROS AND U IS A 
VECTOR STORED UPPER TRIANGULAR MATRIX 


phi ( irphi » icphi ) 

MAXPHI 

IRPHI 

ICPHI 

U(N*(N+D/2) 

N 

PHIUdRPHIiN) 

MPHIU 


INPUT RECTANGULAR MATRIX t TRPHT .LE. MAXPHI 

ROW dimension of phi 

NO. ROWS OF PHI 
NO. COLS OF PHI 

UPPER TRIANGULAR VECTOR STORED MATRIX 
DIMENSION OF U MATRIX (ICPHI. LE-N) 

OUTPUT * RESULT OF PHl*Ur PHIU CAN 

OVERWRITE PHI 

ROW DIMENSION OF PHlU 


COGNIZANT PERSONS? G. J.BIERMAN/M. W»NEAD (JPL» FEB. 19781 


IMPLICIT DOUBLE PRECISION (A-H»0-Z> 
'DIMENSION PHI (MAXPHI * 1 ) *U(1) »PHIU(MPHIU#1> 
DOUBLE PRECISION SUM 

•DO 10 1=1 » IRPHI 
10 PHIU ( I » 1 ) =PHI ( I . 1 ) 

• ,NP2=N+2 

KJS=N*(N+l)/2 
DO '40 L=2»N 
‘ * • J=NP2-L 

iKJS=KJS-vJ 

,UM1=J-1 

- . DO 30 1=1# IRPHI 

- SUM=PHl ( I # 

• IF (J.LE# ICPHI) GO TO 15 
SUM=0.oO 

- JMl'=ICpHI 

15 DO - 20 k=1»»JM1 

20 SUM=SUm+PHI(I»K)*U(KJS+K) 

30 PHIU(I#J)=SUM 
■40 -CONTINUE 

RETURN 

END 


PHIU0040 

PHIU0050 

PHIU0060 

PHIU0070 

PHIU0080 

PHIU0090 

PHIU0100 

PHIU0110 

PHIU0120 

PHIU0130 

PHIU0140 

PHIU0150 

PHIU0160 

PHIU0170 

PHIUQ180 

PHIU0190 

PHIUOPOO 

PHIU0210 

PHIU0220 

PHIU0230 

PHIU0240 

PHIU0250 

PHIU0260 

PHIU0270 

PHIU0280 

PHIU029O 

PHIU0300 

PHIU0310 
PHIU0320 
PHIU0330 
PHIU0340 
PHIU0350 
PHIU0360 
PHIU0370 
PHIU03S0 
PHI UO 390 
PHIU0400 
PHIU0410 
PHIU0420 
'PHIU0430 
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SUBROUTINE RA <R»N»A*MAXA*IA' JA»RA*MAXRA»NRA) 

TO COMPUTE RA=r*A 

WHERE R IS UPPER TRIANGULAR VECTOR SUBSCRIPTED AND OF DIMENSION N 
A HAS JA COLUMNS AND IA ROWS* IF IA.LT.JA THEN THE BOTTOM JA-Ia 
ROWS OF A ARE ASSUMED TO BE IMPLICITLY DEFINED AS THE 
BOTTOM JA-IA ROWS OF THE JA DIMENSION IDENTITY MATRIX. 

ONLY NR A ROWS OF THE PRODUCT R*A ARE COMPUTED. 


& R<N*(N+l>/2) 
N 

A(IA*Ja> 

MAXA 

IA 

JA 

RA(NRA.N) 

maxra 

NRA 


UPPER TRIANGULAR VECTOR STORED INPUT MATRIX 

DIMENSION OF R 

INPUT RECTANGULAR MATRIX 

ROW DIMENSION OF A 

NUMBER OF ROWS IN THE A MATRIX (IA.LE.MAXA) 

NUMBER OF COLUMNS IN THE A MATRIX 
OUTPUT RESULTING RECTANGULAR MATRIX * 

RA=A IS ALLOWED 
ROW DIMENSION OF RA 

NUMBER OF ROWS OF THE PRODUCT R*A THAT ARE COMPUTED 
(NRA. LF. MAXRA) 


COGNIZANT persons: G* J.BIERMAN/M.W.nFAD 

IMPLICIT DOUBLE PRECISION <A-H*0-Z) 
DIMENSION R(1).A(MAXA>I) >RA(MAXRAfl) 
DOUBLE PRECISION SUM 


(JPL» FEB. 1978) 


IJ=IA*(IA+l)/2 


0 IJ=JJUA) 

0 TO RE REMOVED IF JJ(I) IS USED 


10 

15 

20 

30 


DO 30 J=1»JA 
11=0 
DO 2() 1=1 *NRA 

II = II + I 6> II=(I»I)=JJ(I) 

IT IS MORE EFFICIENT TO USE A PRESTOREO VECTOR OF DIAGONALS 
WITH JJ(I)=I*(I+l)/2» AND TO S^T II=JJ(I) AND IJ=JJ<J) 

SUM=0.00 

IF (I.GT.IA) GO TO 15 
IK=II 

DO 10 K=I»IA 

SUM=SUM+R(IK)*A(K»J) 

IK^IK+K 

IF (J.GT.IA.AND.I.LE.J) SUM=SUM+R(IJ+T) 


RA ( I » J) =SUM 
IF (J.GT.IA) IJ=IJ+J 

RETURN 

END 


13 IJ=JJ(J) 


RA000010 
RA000020 
RA 000030 
RA000040 
f RA000050 
RA000060 
RA000070 
RA000080 
RA000090 
RA000100 
RA000110 
RA000120 
RA000130 
RA000140 
RA000150 
RA000160 
RA000170 
RA000I 80 
RA000190 
RA000200 
RA000210 
RA000220 
RA000530 
RA000240 
RA000250 
RA000260 
RA000270 
RA000280 
RA000290 
RA000300 
RA000310 
RAOO032O 
RA000330 
RA000340 
RA000350 
RA000360 
RA000370 
RA000380 
RA000390 
RA000400 
RA000410 
RA000420 
RA000430 
RA000440 
RA000450 
RA000460 
RA000470 
RA000480 
RA000490 
..RA000500 
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SUBROUTINE RANK1 <UIN*U0UT*N*C»V) 

STABLE U-D FACTOR RANK 1 UPDATE 

(U0UT)*D0UT*<U0UT)**T=<UIN)*DIN*<UIN>**T+C*V*V**T 

UIN(N*(N+l)/2) input vector stored positive semi-definite u»d 

ARRAYr WITH D ELEMENTS 5T0RE0 ON THE DIAGONAL 
UOUT(N*(N+l)/2) OUTPUT VECTOR STORED POSITIVE (POSSIBLY) SEMJ- 
DEFINITE U-D RESULT. U0UT=UIN IS PERMITTED 
N MATRIX DIMENSION* N.GE.2 

C INPUT SCALAR. SHOULD BE NON-NESATIVF 

C IS DESTROYED DUPING THE PROCESS 
V(N) INPUT VECTOR FOR RANK ONE MODIFICATION. 

V IS DESTROYED DURING THE PROCESS 

COGNIZANT PERSONS: G. J.BIERMAN/M.W.nEAD (JPL* SEPT. 1977) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION UIN { 1 ) * UOUTd)* V(l) 

DOUBLE PRECISION ALPHA* BETA* S* D* Eps* TST 

DATA EpS/0. DO/* TST/.0625D0/ 

IN SINGLE PRECISION EPSILON IS MACHINE ACCURACY 

TST=l/l6 IS USED FOR RANKl ALGORITHM SWITCHING 

Z=0.D0 

JJ=N*(n+1)/2 
IF (C.GT.Z) GO TO 4 
DO 1 J=1*JJ 

1 UOUT(J)=UIN(J) 

RETURN 

4 NP2=N+2 

DO 70 L=2*N 
J=NP2-L 
S=V(J) 

BETA=C*5 

D=UlN<<JJ>+BETA*S 
IF (D.GT.EPS) GO TO 30 
IF (o.GE.Z) GO TO 10 

5 WRITE (6*100) 

RETURN 

10 JU— *J»J— J 

WRITE (6*110) 

DO 2o K=1*J 

20 UOUT ( JJ+KJ=Z 

GO To 70 

30 BETA=BETA/D 

ALPHa=UIN(JJ)/D 
c=alpha*c 
UOUT ( JJ) -D 


RANK1010 
RANK1020 
RANK1030 
RANK1040 
RANK1050 
RANK1060 
RANK1070 
RANK1080 
RANK1090 
RANK1100 
RANK1110 
RANK1120 
RANK1130 
RANK1140 
RANK1150 
RANK 1160 
RANKl 170 
RANK1180 
RANKl 190 
RANK1200 
RANK1210 
RANK1220 
RANK1230 
RANK1240 
RANK1250 
RANK1260 
RANK 1270 
RANK1280 
RANK1290 
RANK1300 
RANK1310 
RANK1320 
RANK1330 
RANK1340 
RANK1350 
RANK1360 
RANK1370 
RANK1380 
RANK1390 
RANK1400 
RANK1410 
RANK1420 
RANK1430 
RANK1440 
RANK1450 
RANK1460 
RANK1470 
RANK1480 
RANK1490 
RANK1500 
RANK1510 
RANK1520 
RANK1530 
RANK1540 
RANK1550 
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pagb ® 


IF ( ALPHA. LT.TST) GO TO 50 
DO 40 I=1*JM1 

V ( I ) =V c I ) -S*UIN ( J J+I > 

UOUT < -J J+I ) =BETA*V ( I ) +UIN ( J J+I ) 

GO TO 70 
DO 60 I~1 » JM1 

0=Vd)“S+UIN( JJ+I) 

UOUT ( JJ+I )=ALPHA*UIN ( JJ+I ) +BETA+V ( I ) 
V<I>=D 
CONTINUE 

UOUT ( 1 ) =UIN ( 1 ) +C*V ( 1 ) **2 
RETURN 


100 FORMAT (IHOflOX*** * * ERROR RETURN DUE TO 
1PUTED DIAGONAL IN RaNKI * * **> 

110 FORMAT (1H0»10X»»* * * note: u-d RESULT is 
end 


40 

50 


60 

70 


RANK1560 
RANK1570 
RANK1580 
RANK1590 
RANK 1600 
RANK1610 
RANK 1620 
RANK1630 
RANK1640 
RANK 1650 
FANK1660 
RANK 1670 
RANK1680 
RANK1690 

A COMPUTED NEGATIVE COMR ANK1700 

RANK1710 

SINGULAR * * *') RANK1720 

RANK1730 
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SUBROUTINE PCOLRD (S *MAXS* IRS * JCS*NPSTRT»NP *FM* RW*ZW* V * SGST AR) 

TO ADD IN PROCESS NOISE EFFECTS INTO THE SOUARE ROOT 

information filter* and to generate weighting coefficients 

FOR SMOOTHING. IT IS ASSUMED THAT VARIABLES X(NPSTRT)* 
X(NPSTRT+i)*...*X(NPSTRT+NP“l) ARE COLORED NOISE AND THAT 
EACH COMPONENT SATISFIES A MODEL EQUATION OF THE FORM 
X(SUB) (J+1)=EM*X<SUB) (J)+W(SUR) (J) . FOR DETAIL5* SEE 
♦FACTORIZATION methods FOR DISCRETE sequential estimation** 
G.J.BIERMAN* ACADEMIC PRESS (1977) 

FOR SMOOTHING* REMOVE THE COMMENT STATEMENTS ON THE 3 LINES 
OF ’SMOOTHING ONLY* CODE* THE SIGNIFICANCE OF THE SMOOTHING 
MATRIX IS EXPLAINED IN THE FUNCTIONAL DESCRIPTION. 


S ( IRS * JCS) 


MAXS 


IRS 

JCS 


NPSTRT 

NP 

EM(NP) 

RW(NP) 

ZW(NP) 


V ( IRS) 
SGSTAR(NP) 


INPUT SQUARE ROOT INFORMATION ARRAY. OUTPUT COLORED 

noise array housed here too. if there is smoothing* 

NR ADDITIONAL ROWS ARE INCLUDED IN S 
ROW DIMENSION OF S. IF THERE ARE SMOOTHING COMPUTA- 
TIONS IT IS NECESSARY THAT MAXS .GE. IRS+NP BECAUSE 
THE BOTTOM NP ROWS OF S HOUSE THE SMOOTHING 
INFORMATION 

NUMBER OF ROWS OF S (.LE* NUMBER OF FILTER VARIABLES) 
(IRS.GE.2) 

number of columns of s (Equals number of filter 

VARIABLES + POSSIBLY A RIGHT SIDE)* WHICH CONTAINS 
THE DATA EQUATION NORMALIZED ESTIMATE (JCS.GE.l) 
LOCATION OF THE FIRST COLORFD NOISE VARIABLE 
(l.LE. NPSTRT. LE. JCS) 

NUMBER of CONTIGUOUS COLORED NOISE VARIABLES 
COLORED NOISE MAPPING COEFFICIENTS 
(OF EXPONENTIAL FORM* EM=EXP (-DT/TAIJ) ) 

RECIPROCAL PROCESS NOISE STANDARD DEVIATIONS 
(MUST BE POSITIVE) 

ZW=RW*W-ESTIMATE (PROCESS NOISE ESTIMATES ARE 
GENERALLY ZERO MEAN). WH^N 7W=0 ONE CAN OMIT THE 
RIGHT HAND SIDE COLUMN. 

WORK VECTOR 

VECTOR OF SMOOTHING COEFFICIENTS. WHEN THE SMOOTHING 
CODE IS COMMENTED OUT SGSTAR IS NOT USED. 


COGNIZANT PERSONS: G. J.BIERMAN/M.W.nEAD (JPL* FEB. 1978) 
IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION S(MAXS*JCS)*EM(NP)*RW(NP)*ZW(NP)» V( IRS) *SGSTAR ( 1 > 
DOUBLE PRECISION ALpHA*SlGMA*BETA*GAMMA 


ZERO=O.DO 

ONE=1.DO 

npcol=npstrt 


0 COL NO OF COLORED NOISE TERM TO BE OPERATED ON 


DO 70 JC0LRD=1*NP 

ALPHA=-RW ( JCOLRD ) *EM ( JCOLRP) 
sigma=alpha**? 

DO lo K=1*IRS 

v(k)=s(k*npcol) (3 FIRST IPS 


ELEMENTS of householder 

104 


RCOLR010 
RCOLR020 
RCOLR030 
RCOLR040 
RCOLR050 
RCOLR060 
RCOLR070 
RCOLR080 
RCOLR090 
RCOLR100 
RCOLR110 
RC0LR12Q 
RCOLR130 
RC0LR140 
RCOLR150 
RC0LR160 
RCOLR170 
RC0LR180 
RCOLR190 
RCOLR200 
RCOLR210 
RCOLR220 
RCOLR230 
RCOLR240 
RCOLR250 
RCOLR260 
RC0LR270 
RCOLR280 
(NP.GE.DRCOLR290 
RC0LR300 
RC0LR310 
RCOLR320 
RC0LR330 
RC0LR340 
RCOLR350 
RCOLR360 
RC0LR370 
RCOLR380 
RC0LR390 
RCOLR400 
RCOLR410 
RCOLR420 
RCOLR430 
RCOLR440 
RCOLR450 
RCOLR460 
RC0LR470 
RCOLR480 
RC0LR490 
RCOLR500 
RCOLR510 
RC0LRS20 
RC0LR530 
RCOLR‘540 
RCOLRB50 
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C 

10 


* * 
* * 


* * 
* * 
20 
30 


40 


50 
* * 

* * 
60 
* * 


* * 
70 


TRANSFORMATION VECTOR 

SIGMA=SIGMA+V (K ) **2 
SIGMA=DSQRT(SIGMA) 

ALPHA=ALPHA-SIGMA 13 EAST ELEMENT OF HOUSEHOLDER 

transformation VECTOR 

* * * * 

SGsTAR(JCOLRD)=sIGMA 0 USED for smoothing only 
* * * * 

BETA=0NE/(SIGMA*ALPHA) 13 householders +beta*v*v**t 

HOUSEHOLDER TRANSFORMATION DEFINED , Now APPLY IT TO S> I.E .60 
DO 6o KOL=l»JCS 

IF (K0L.NE.NPCOL) GO TO 30 
GAmMA= RW(JCOLRD)*AUPHA*RETA 
* * * * 

S(IRS+JCOLRD»NPCOL)=RW(JCOLRD)+GAMMA*ALPHA q smoothing 
* * * * 

DO 20 K=1»IRS 

S (K * NPCOL) = GAMMA*V(K ) 

GO TO 60 
GAMMA=2ERO 

if <kol.eq*jcs) GAMMA=2W ( JCOLRD ) * alpha 
IF ZW ALWAYS ZERO* COMMENT OUT the ABOVE IF TEST 


RC0LR560 
RCOLR570 
RCOLR530 
RCOLR590 
RCOLR600 
RCOLR610 
RCOLR620 
RCOLR630 
RCOLR640 
LOOPRCOLR650 
RC0LR660 
RCOLR670 
RCOLR680 
RCOLR690 
ONLY RCOLR7O0 
RCOLR710 
RCOLR720 
RC0LR730 
RC0LR740 
RC0LR750 
RC0LR760 
RC0LR770 
RCOLR780 


DO 40 KS*IRS 

GAMMA=GAMMA+S{KtKOL)*V<K) 

gamma= gamma*beta 

DO 50 K=1»IRS 

S(K*K0U=5<K*K0L)+GAMMA*VtK) 

* * * * 

s< irs+jcolRd*kol)=gamma*alpha 
* * * * 

CONTINUE 
* * * * 


RC0LR790 

RCOLR800 

RC0LR810 

RC0LR820 

RC0LR830 

RCOLR840 

RC0LR850 

0 FOP SMOOTHING ONLY RC0LR860 

RCOLRB70 
RC0LR880 
RC0LR890 
RCOLR900 

RC0LR910 
RCOLR920 
RC0LR930 
RC0LR940 
RC0LR950 
RCOLR960 
RC0LR970 


5 < I RS+JCOLRO * JCS > =S ( IRS+ JC OLRD * jCS ) +ZW ( UCOLRD ) 

the above is for smoothing ONLY 

if ZW is always ZERO, comment out the above statement 
* *' * * 

NPCOL=NPCOL+l 

RETURN 

END 
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* SUBROUTINE RINCON ( Rl N , N » ROUT;, CNB ) 

TO COMPUTE THE INVERSE OF THE UPPER TRIANGULAR VECTOR STORED 
INPUT MATRIX RIN-AND STORE THE RESULT IN ROUT. (RIN=ROUT IS 
PERMITTED) AND TO COMPUTE A- CONDITION NUMBER ESTIMATE. 
CNB=FR0B.N0RM(R)*FR0B.N0RM(R**-1) , 

The FRoRENIUS norm is the square root of the sum of squares 
of THE ELEMENTS. THIS CONDITION NUMBER BOUND IS USED AS 
AN UPPER BOUND AND IT ACTS AS, A LOWER BOUND ON THE ACTUAL 
CONDITION NUMBER OF THE PROBLEM. (SEE THE BOOK ‘SOLVING LEAST 
SQUARES’* BY LAWSON AND HANSON) 

IF RIN IS SINGULAR, RINCON COMPUTES THE INVERSE TO THE LEFT OF 
THE FIRST ZERO DIAGONAL. A MESSAGE 1^ PRINTED AMD THE CONDITION 
NUMBER BOUND COMPUTATION IS ABORTED. 

RIN(N*(N+l)/2) INPUT VECTOR STORED UppfR TPI ANGULAR MATRIX 
N DIMENSION OF R MATRICES, N.6E.2 

ROUT<N*(N+l)/2) OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX 

inverse (rin=rout is permitted) 

CNB CONDITION NUMBER BOUND. IF C IS THE CONDITION 

NUMBER OF RIN, THEN CMB/N.LF.CiLE.CNB 

COGNIZANT PERSONS: G* J.B TERM AN/M, w«f|E AD ( JPL» FEB . 1978) 


IMPLICIT DOUBLE PRECISION <A-H»0-Z) 

DOUBLE PRECISION RNM» DINV» SUM»RNMOUT 
DIMENSION RIN(l), ROUT(l) 

Z=0.D0 
ONE=1.DO 
NT0T=N*<N+l>/2 

RNM=Z 

DO 10 J=l»NTOT 
10 RNM=RNM+RIN<J)**2 

REPLACE CALL UTlNV (RIN, N, ROUT) BY UTINV CODE 

IF (RlNd) .NE.Z) GO TO 20 
J=1 

WRITE (6,100) J,J 
RETURN 

20 ROUT(l)=ONE/RIN(l) 

C 

JJ=1 

DO 50 J=2,N 
J J0Ld = <J*J 
J 

IF (RlN(JJ) .NE.Z) GO TO 30 
WRITE (6,100) J»J 
RETURN 

C 106 


RTMC0010 
RINC0020 
RINC0030 
RINC0040 
RTNC0050 
RINC0060 
RIMC0070 
RINCOOSO 
RTNC0090 
RIMC0100 
RINC0110 
RINC0120 
RTNC0130 
RINC0140 
RINC0150 
RINC0160 
RTNC0170 
RTNC0180 
RINC0190 
RTNC0200 
RINC02in 
RINC0220 
RTNC0230 
RINC0240 
RINC0250 
RINC0260 
RINC0270 
RINC0280 
RINC0290 
RINC0300 
RTNC0310 
RINC0320 
RINC0330 
RTNC0340 
RTNC0350 
RINC0360 
RTNC0370 
RTNC0380 
RINC0390 
RTNC0400 
RINC0410 
RINC0420 
RINC0430 
RTNC0440 
RTMC0450 
RINC0460 
RINC0470 
RINC0480 
RTNC0490 
RINC0500 
RINC0510 
RINC0520 
RTNC0530 
R INC 0540 
RINC0550 
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C 

C 

C 


30 DINV=ONE/RIN(JJ) 

ROUT { JJ)=DINV 
11=0 
IK=1 
JM1=J-1 
DO 5o 1=1 » JM1 
II=II+I 
IKxII 
SUM=Z 

DO 40 K=IfJMl 

SUM=SUM+ROUT ( I K ) *R IN ( J JOLD+K ) 
40 IK=IK+K 

50 ROuT(JJOLD+I)=-SUM*DINV 


RNMOUTxZ 
DO 60 J=l*NTOT 

60 RNM0UT=RNM0UT+R0UT(J)**2 

RNM=DSqRT(RNM*RNMOUT> 

CNB=RNM 

WRITE (6» 110) RNM 
RETURN 


RINC0560 

RINC0570 

RTMC0580 

RINC0590 

RINC0600 

RINC0610 

RTNC0620 

RINC0630 

RINC0640 

RINC0650 

RTNC0660 

RTNC0670 

RTNC0680 

RTNCO690 

RTNC0700 

RINC0710 

RTNC0720 

RTNC0730 

RINC0740 

R1NC0750 

RINC0760 

RTNC0770 

RINC0780 

RINC0790 

RINC0800 

RINC0810 


100 FORMAT (lH0»10Xr '* * * MATRIX INVERSE COMPUTED ONLY UP TO BUT NOT RTNCO02O 
1INCLUDING COLUMN* »I4»* * * * MATRIX DIAGONAL * *I4»* IS ZERO * * **RTNC0830 
2) RINC0840 

110 FORMAT ( 1H0 » 5X » ’CONDITION NUMBER ROUMD=* pDl8.lO»2X t » CNB/N«LE.C0N0ITRINC0850 
II ON NUMBER»LE.CNB* f/) RINC0860 

End RTNC0870 
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SUBROUTINE RI2C0V (rINV*N»SIG»COVOUT»KROW»KCOL> RT2C0010 

RT2COC120 

TO COMPUTE THE COVARIANCE MATRIX AND/OR THE STANDARD DEVI ATIONSRI2C0030 


C 

c 


OF A VECTOR STORED UPPER TRIANGULAR SQUARE ROOT COVARIANCE 
MATRIX. the output covariance matrix IS ALSO vector STOREO. 


RINV<N*(N+l)/2) 


N 

SIG(N) 

C0V0UT(N*(N+l)/2) 
KROW .GT.O 


.LT.O 


.EQ.O 


KCOL 


INPUT VFCTOR STORED UPPER TRIANGULAR 
COVARIANCE SQUARE ROOT. (RINV=RINVERSE 
IS THE INVERSE OF THE SRIF MATRIX) 
dimension OF the rinv matrix. N.GE.2 
output vector of standard deviations 
output vector STORED COVARTAnCF MATRIX 
(COVOUT = RINV is ALLOWED) 

COMPUTES THE CoVARTANCE AND SIGMAS 
CORRESPONDING TO THE FIRST KROW VARIABLES 
OF THE RINV MATRIX. 

COMPUTES ONLY ^E SIGMAS OF THE FIRST KROW 
VARIABLES OF THF RINV MATRIX. 

RINV. 

MO COVARIANCF. PUT ALL SIGMAS (F»G, USE 
N ROWS OF RINV)* 

NO* OF COLUMNS OF COVONT THAT ARE COMPUTED 
IF KCOL*LE.O THEN KCOL=KROW. IF KROW.LE.O 
THIS INPUT IS IGNORED. 


COGNIZANT PERSONS: G. J..RIERMAN/M* W.NEAD <JPL* MARCH 1978) 

IMPLICIT DOUBLE PRECISION (A-H.O-2) - 

DOUBLE PRECISION SUM 

DIMENSION RINV(l). SIG(1) . cOVOUT(l) 

ZERO=O.DO 

LIM=N 

KK0L=KC0L 

IF (KKOL.LE.O) KKOL=KROW 
IF CKROW*NE.O) LIM=IABSCKR0W) 

*** COMPUTE SIGMAS 

IKS=0 

DO 2 Uzl.LIM 
IKS=IKS+U 
SUM=ZERO 
IK=IkS 
DO 1 K=U»N 

SUm=SUM+RINV<IK)**2 

IKrlK+K 

SIG<J)=DSQRT(SUM) 

IF (KROW.LE.O) RETURN 

*** COMPUTE COVARIANCE 

JJ=0 

NM1=LIM 

IF (KROW.EQ.N) NM1=n- 1 
DO 10 J=1»NM1 
JJ=Jj+J 

COVOUT (JJ)=SIG<J)**2 


RT2C0040 

RI2C0050 

RT2C0060 

RI2C0070 

RI2COOSO 

RI2C0090 

RT2C0100 

RI2C0110 

RT2CO120 

RT2C0130 

RI2C0140 

RT2CO150 

RI2C0160 

RT2C0170 

RI2CO1Q0 

RT2C0190 

RI2C0200 

RI2CO210 

RI2CO220 

RI2C0230 

RI2CO240 

RI2C0250 

RI2C0260 

RI2CO270 

RI2C0280 

RI2C0290 

RI2C0300 

RT2CO310 

RI2C0320 

RI2C0330 

RT2C0340 

RI2C0350 

RI2C0360 

RI2C0370 

RT2CO380 

RI2C0390 

RI2C0400 

RI2C0410 

RI2C0420 

RI2C0430 

RI2C0440 

RI2C0450 

RI2C0460 

RT2C0470 

RI2CO480 

RI2C0490 

RI2C0500 

RI2CO510 

RI2C0520 

RT2C0830 

RI2C0540 

RI2C0550 
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IJS=JJ+J 

jpi=j+i 

DO lo I=JPl»KKOL 
IK=IJS 
IMJ=I-J 

sum=zero 

DO 5 K=I»N 
IJK=IK+IMJ 

SUM=5UM+R INV(IK)*RINV(IJK) 

5 IK=IK+K 

COVOUT(IJS)=SUM 
10 IJS=IJS+I 

IF (KROW.EQ.N) COVOUT(JJ+N)=SIG(N)**2 

RETURN 

END 


RT2C0560 

RI2C0570 

RI2CO580 

RT2CO590 

RT2C0600 

RI2CO610 

RI2C0620 

RI2C0630 

RI2CO640 

RI2C0650 

RI2CO660 

RI2C0670 

RI2CO680 

RI2CO690 

RI2C0700 

RT2CO710 
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SUBROUTINE R2 A ( R » LR » NAMR » A » I A » LA » MAMA ) 


TO PLACE THE TRIANGULAR VECTOR STORED MATRIX R INTO THE 
MATRIX A AND TO ARRANGE THE COLUMNS TO MATCH THE DESIRED 
NAMa PARAMETER LIST. NAMES IN THE NAMA LIST THAT DO NOT 
CORRESPOND TO ANY NAME IN NAMR HA^E 7ER0 ENTRIES IN THE 
CORRESPONDING A COLUMN. 


R(Lr*(LR+1)/2) 

LR 

NAMR(L) 

A(LR»LA) 

IA 

LA 

>NAMa<LA) 


INPUT UPPER TRIANGULAR VECTOR STOREO ARRAY 
OtMENSlON OF R 

PARAMETER NAMES ASSOCIATED WITH R 
MATRIX TO house the rearranged R MATRIX 
ROW DIMENSION OF A, IA.GE.LP 

mo. of parameter names associated with the 

OUTPUT A MATRIX 

PARAMETER names for the output a matrix 


COGNIZANT persons: G.J.BlERMAN/MtW.NEAD (JPL» SEPT. 1976) 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
DIMENSION R<1) 'NAMR(l) 'AQA'l ) * NA M A ( 1 ) 

ZERO = 0 • 

DO 5 J=1»LA 
DO 5 K=lfLR 

5 A(K.J)=ZERO Q ZERO AtLR.LA) 

DO 40 J=1»LA 
DO lo I = 1 * LR 

IF (NAMR(I) .EQ.nAMA(J) ) GO TO 20 
10 CONTINUE 

GO TO 40 
20 JJ=I*(l-l>/2 
DO 3o K=1»I 
30 A(K» J)=R(JJ+K> 

40 CONTINUE 

RETURN 

END 


R2AOOO10 
R2A00020 
, R2A00030 
R2A00040 
R2A00050 
R2A00060 
R2A00070 
R2A00080 
R2A00090 
R2A00100 
R2AOOU0 
R2A00120 
R2A00130 
R2A00140 
R2A00150 
R2A00160 
R2A00170 
R2A00180 
R2A00190 
R2A00200 
R2A00210 
R2A00220 
R2A00230 
R2A00240 
R2A00250 
R2A00260 
R2A00270 
R2A00280 
R2A00290 
R2A00300 
R2A00310 
R2A00320 
R2A00330 
R2A00340 
R2A00350 
R2A00360 
R2A00370 
R2A00380 


no 



no o o o oooooonoonnnooooooooooooooooonnooooo 


ORIGINAL RAGE IS 
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1 

3 

100 


SUBROUTINE R2RA (R*nR>NAM*RA'NPA»NAMA) 

TO COPY THE UPPER LEFT (LOWER RIGHT) PORTION OF A VECTOR 
STORED UPPER TRIANGULAR MATRIX R INTO THE LOWER RIGHT 
(UPPER LEFT) PORTION OF A VECTOR STORED TRIANGULAR 
MATRIX RA. 


R(Nr*(NR+1>/2) 

NR 

NAM (NR) 

RA(NRA*(NRA+1 )/2) 
NR A 


NAMa(NRA) 


INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 

DIMENSION OF R 

NAMES ASSOCIATED WITH R 

THIS INPUT NAMELIST IS DESTROYED 

OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX 

IF NRA=0 ON INPUT* THEN NAMA(l) SHOULD HAVE 

THE FIRST NAME OF THE OUTPUT NAMELIST. 

IN THIS CASE THE NUMBER OF NAMES IN NA M A AND 
NRA WILL RE COMPUTED. THE LOWER FIGHT FLOCK 
OF R WILL BE TH^ UPPER LEFT BLOCK OF RA, 

IF nra=last name OF the UPPER left block 
THAT IS TO BE MOVED » THEN THIS UPPER 
BLOCK IS TO BE MOVED TO THE LOWER RIGHT 
CORNER OF RA. WHEN USED IN THIS MODE NRA=NR 
ON OUTPUT. 

NAMES ASSOCIATED WITH RA 


IF NRA=0 ON INPUT* THEN NAMAd> SHOULD HAVE THE FIRST NAME OF THE 
OUTPUT NAMELIST AND THE NUMBER OF NAMES IN MAMA IS COMPUTED. 

THE LOWER RIGHT BLOCK OF R WILL BE THE UPPER LEFT BLOCK OF RA. 

IF NRA=LAST NAME OF THE UPPER LEFT BLOCK THAT IS TO BE MOVED* 
then The upper block is TO be MOVED To the LOWER RIGHT POSITION. 
WHEN USED IN THIS MODE NRA=nR ON OUTPUT. 

THE NAMES OF THE RELOCATED BLOCK are ALSO MOVED. THE RESULT 
CAN COINCIDE WITH R AND NAMA WITH NAM. 


R2PA0ni0 

R2RA0P2O 

RPRA0030 

R2RA0040 

R2RA0050 

R2RA0060 

R2RAD07R 

RPRAOnBO 

R2RA0090 

R2RA0100 

R2RA0110 

R2RA0120 

R2RA0130 

R2RA0140 

R2RA0150 

R2RA0I60 

R2RA0170 

R2RA018O 

R2PA0190 

R2RA0200 

R2RA0210 

R2PA0220 

R2PA0230 

R2RA0240 

R2RA0250 

R2RA0260 

R2RA0270 

R2RA0280 

R2RA029O 

R2RA0300 

R2RA031O 

R2RA0320 

R2RA0330 

R2RA0340 


COGNIZANT PERSONS? G. J.BIERMAN/M. W.NEAD (JPLf SEPT, 1976) 

IMPLICIT DOUBLE PRECISION (A-H»0-Z> 

DIMENSION R ( 1 ) * RA( 1 ) * NAM<1)* NAMA(l) 

logical IS 

IS=. FALSE. 

LOCN=NAMA(l) 

IS=FALSE CORRESPONDS TO MOVING UPPER LFT . CORNER OF R TO 
LOWER RT. CORNER OF RA 
IF (NRA.EQ.O) GO TO 1 
LOCN=NRA 
IS=.TRUE. 

IS=TRUE CORRESPONDS TO MOVING LOWER LFT. CORNER OF R TO 
UPPER RT. CORNER OF RA 
DO 3 1=1* NR 

IF (NAM ( I > .EG.LOCN) GO TO 4 
CONTINUE 
WRITE (6*100) 

FORMAT (1H0*20X* ’NAMA(l) NOT IN NAMEt-lST OF R MATRIX*) 


R2RA0350 

R2RA0360 

R2RA0370 

R2PA0380 

R2RA0390 

R2RA0400 

R2RA0410 

R2RA0420 

R2RA0430 

R2RA0440 

R2RA0450 

R2RA0460 

R2RA0470 

R2RA0480 

R2PA0490 

R2RA0500 

R2RA0510 

R2RA0520 

R2RA0530 

R2RA0540 

R2RA0S50 
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RETURN 

4 K=I 

KMl=K-l 

IF (IS) GO TO 15 

IJS=K*(K+l)/2-l 

NRA=NR-K+1 

IJA=0 

KOLA=0 

DO 10 KOL=K,NR 
koua=kola+i 

NAM A ( K0L-KM1 ) =N AM ( KOL ) 

DO 5 IR=1»K0LA 

ua=ija+i 

5 RA(IJA)=R(IJS+IP> 

10 IJS=lJS+KOL 

RETURN 

15 IJ=K*CK+l)/2 

IJA=Nr*(NR+1)/2 

L=NR-kM1 

KOL=K 

DO '25 KOLA=NR>L#-l 
IJS=IJA 

NAMA(KOLA)=NAM(KOL) 
DO 20' IR=KOLA»Lr-i 
RA(IJS)=R(IJ) 

20 IJ=IJ-1 

ija=ija~kola- 

.25 KOL=KOL-l 

nra=nr 

RETURN 

END 


R2RA0560 

R2PA0570 

R2RA0580 

R2RA0590 

R2RA0600 

R2RAO610 

R2RA0620 

R2RA0630 

R2RA0640 

R2RA0650 

R2RA0660 

R2RA0670 

R2RA0680 

R2RA0690 

R2RA0700 

R2RA0710 

R2RA0720 

R2RA0730 

R2RA0740 

R2RA0750 

R2RA0760 

R2RA0770 

R2RA0780 

R2RA0790 

R2RA0800 

R2RA0810 

R2RA0820 

R2RA0830 

R2RA0840 

R2RA0850 

R2RA0860, 

R2RA0870 

R2RA0880 

R2RA0890 

R2RA0900 

R2RA0910 
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C 


c 


c 


c 

c 

c 


SUBROUTINE RUDR (RIN.N »ROUT » IS) 


FOR N.GT.O THIS SUBROUTINE TRANSFORMS AN UPPER TRIANGULAR VECTOR 
STORED SRIF MATRIX TO U^D FORM. AND WHEN N.LT.O THE U-D VECTOR 
STORED ARRAY IS TRANSFORMED TO A VECTOR STORED SRIF ARRAY 


RIN( (N+l)*(N+2>/2) 
ROUT ( (n+1>*<N+2)/2) 

N 

N » GT • 0 
N.LT.O 
IS = Q 

IS = 1 


input vector stored srif or u-d array 
output is THE CORRESPONDING u-d or srif 
ARRAY (RIN=ROUT IS PERMITTED) 

ABS(N)= MATRIX DIMENSION .GE.2 
THE (INPUT) SRIF ARRAY IS (OUTPUT) 

IN U-D FORM 

THE (INPUT) U-D ARRAY IS (OUTPUT) 

IN SRIF FORM 

THERE IS NO RT. SIDE OR ESTIMATE STORED In 
COLUMN N+l » AND RIN NEED HAVE ONLY 
N COLUMNS ► I.E. RlN<N*(N+l)/2) 

THERE IS A RT. SIDE INPUT TO THE SRIF AND 
AN ESTIMATE FOR THF U-D ARRAY. THESE RESIDE 
IN COLUMN N+l. 


THIS SUBROUTINE USES SUBROUTINE RINCON 


COGIZANT PERSONS G.J.BIERMAN/M.W.NEAD ( JPL » FER.197B) 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
DIMENSION RIN(1). ROUT(l) 


ONE= 1 . DO 

NP1= IS + IABS(N) 

JJ=1 (3 initialize DIAGONAL index 

IDIMR= NP1*(NP1 +l)/2 
IF (IS.EQ.O) GO TO 5 
RNN=RIN<IDIMR) 

RIN(IDlMR)=-ONE 

5 IF (N.LT.O) GO TO 30 

CALL RlNCON(RlN^NPlfROUTrCNR) 

ROUT ( 1 ) = ROUT ( 1 ) **2 
DO 20 J=2.N 
S-ONE/ROUT ( JJ+J) 

R0UT(JJ+J)= ROUT<JJ+J)**2 

JM1=J-1 

DO 10 

10 ROUT ( JJ+I ) = ROUT ( JJ+I ) *S 

20 JJ=JJ+ J 
GO TO 70 


30 NN=-N R NN=NEGATIVE N 

ROUT ( 1 ) = SORT (RIN ( 1 ) ) 

*** SOME MACHINES REOUIRE DSORT FOR DOURLE PRECISION 
00 50 J=2rNN 

R0UT(JJ+J)= SQRT(RIN(JJ+J) ) 


RUDR0010 

RUDR0020' 

RUDR0030 

RUDR004C 

RUDR0050 

RUDR0060 

RUDR0070 

RUDR0080 

RUDR0090 

RUDR0100 

RUDROliO 

RUDR0120 

RUDR0130 

RUDR0140 

RUDR0150 

RUDR0160 

RUDR0170 

RUDR0180 

RUDR0190 

RUDR0200 

RUDR0210 

RUDR0220 

RUDR0230 

RUDR0240 

RUDR0250 

RUDR0260 

RUDR0270 

RUDR0280 

RUDR0290 

RUDR0300 

RUDR0310 

RUDR0320 

RUDR0330 

RUDR0340 

RUDR0350 

RUDR0360 

RUDR0370 

RUDR0380 

RUDR0390 

RUDR0400 

RUDR0410 

RUDR0420 
RUDR0430 
RUDR0440 
RUDR0450 
RUDR0460 
RUDR 0470 
RUDR0480 
RUDR0490 
RUDR0500 
RUDR0510 
RUDR0520 
RUDR0530 
RUDR0540 
RUDR0550 
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S=ROUT(JJ+J) 

RIJDR0560 


JM1=J-1 

RUDR0570 


DO 40 I=1»JM1 

RUDR0580 

40 

R0UT(JJ+I)= RIN(«JJ+I ) *5 

RUDR0590 

50 

vJJ=JJ + J 

RUDR0600 

60 

CALL RINCON (ROUT »NPl fROUT »CMB) 

RUDR0610 



RUDR0620 

70 

IF (IS.EQ.l) RIN(IDIMR)=RNN 

RUDR0630 


RETURN 

RUDR0640 


END 

RUDR0650 
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SUBROUTINE SFU ( FEL * IROW * JCOL ♦ NF * U * N » Fu * MAXFlJ * IFU * JDI AG ) 

TO COMPUTE FU(IfU*N)=F*U WHERE F IS SPARSE AND ONLY THE 
NOn-ZERO ELEMENTS ARE DEFINED AND U IS VECTOR STORED* 

UPPER TRIANGULAR WITH IMPLICITLY DEFINED UNIT DIAGONAL 
ELEMENTS 

FEL(NF) VALUES OF THE MON-ZERO ELEMENTS OF THE F MATRIX 

irow (nf > row indices of the f elements 

JCOL(Nf) column indices of the F ELEMENTS 
F(IR0W(K) * JCOL ( K ) ) =FEL ( K ) 

NF NUMBER OF NON-ZERO ELFMENTS OF the F MATRIX 

U(N*(N+l)/2) UPPER TRIANGULAR* VECTOR STORED MATRIX WITH 
IMPLICITLY DEFINED UNIT DIAGONAL ELEMENTS 
(U(J*J) ARE NOT* IN FACT* UNITY) 

N DIMENSION OF U MATRIX 

FU(IFU,N> OUTPUT RESULT 

MAXFU ROW DIMENSION OF FU MATRIX 

IFU NUMBER OF ROWS IN FU* 

( IFU. LF* MAXFU. AND. IFU. OE. MAX (IROW(K) ) * K=l* .. . *NF* 
I.E. FU MUST HAVE AT LEAST AS MANY ROWS AS DOFs F. 
ADDITIONAL ROWS OF FU COULD CORRESPOND To ZERO 
ROWS OF F. 

jdiag(n) diagonal element indices of a vector stored 

UPPER TRIANGULAR MATRIX, 

I.E, JO I AG ( K ) =K* ( K+ 1 ) /2= JDI AG ( K-l ) +K 

COGNIZANT persons: G.J.BIERMAN/M.W.NEAD (JPL* FEB.1978) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION FEL(NF) *U< 1) »FU(MAXFU»N) » IROW(NF) »JCOL(NF> *JDIAG(N) 
ZERO=0,D0 

* * * * INITIALIZE FU 
DO 10 J=1*N 
DO lo 1=1* IFU 
10 FU(I*J)=ZER0 

IF MAXFU-IFU* IT IS MORE EFFICIENT TO REPLACE THIS LOOP BY 

DO 10 IJ=1*1FUN Q IFUM— IFU*N 

10 FU(IJ*l)=ZERO 

DO 30 NEL=1*NF 

NEL REPRESENTS THE ELEMENT NUMBER OP THE F MATRIX 
I=IR0W(NEL) 

J=JCOL(NEL) 

FIJ=FEL(NEU 

FU(I,J)=FU(I,J)+FIJ 

THIS ACCOUNTS FOR THE IMPLICIT UNIT DIAGONAL U MATRIX 
ELEMENTS. WHEN NON-UNIT DIAGONALS ARE USED* DELETE 
THE ABOVE LINE AND USE J INSTEAD OF JPl RFLOW 

IF (J.EQ.N) GO TO 30 

WHEN IT IS KNOWN THAT THE LAST COLUMN OF F IS ZERO 
THIS * IF ’ TEST MAY BE OMITTED 
JP1=J+1 


SFU00010 
SFU00020 
SFU00030 
SFIJ00040 
SFU00050 
SFU00O60 
SFU00070 
SFU00080 
SFU00090 
SFU00100 
SFU00110 
SFU001 20 
SFU00130 
sFuoom-o 
SFU00150 
SFU00160 
SFU00170 
SFU00180 
SPU00190 
SFUOOPOO 
SFU00210 
SFU00220 
SFU00230 
SFU00240 
SFU00250 
SFU00260 
SFU00270 
SFU00280 
SFU00290 
SFU00300 
SFU00310 
SFU00320 
SFU00330 
SFU00340 
SFU00350 
SFU00360 
SFU00370 
SPU00380 
SFU00390 
SFU00400 
SFU00410 
SFU00420 
SFU00430 
SFU00440 
SFU00450 
SFU00460 
SFU00470 
SFU00480 
SFU00490 
SFU00500 
SPU00510 
SFU00520 
SPU00530 
SFU0OS40 
SFUOOS50 
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IK=JdIAG ( J) +J SFU00F60 

DO 2o K=JP1»N SRJ00570 

FU(I»K)=FU(IrK)+FIJ*U(lK) SFU00580 

20 IK=IK+K SFU00590 

30 CONTINUE SF1I00600 

c SFIJ00610 

RETURN SFU00820 

END SRJ00630 
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SUBROUTINE TDHHT ( S * MAXS * IRS * JCS * JSTART » JSTOP * V ) 

tdhht transforms a RECTANGULAR DOUBLE subscripted matrix s 

TO AN UPPER TRIANGULAR OR PARTIALLY UPPER TRIANGULAR FORM 
BY THE APPLICATION OF HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS. 
IT IS ASSUMED THAT THE FIRST * JSTART*-! COLUMNS OF S ARE 
ALREADY TRIANGULARIZED. THE ALGORITHM is described IN 

•factorization methods for discrete sequential estimation* 

BY G.J.BIERMAN* ACADEMIC PRESS* 1977 

S( IRS* JCS) INPUT (POSSIBLY PARTIALLY) TRIANGULAR MATRIX. THp 
OUTPUT (POSSIBLY PARTIALLY) TRIANGULAR RESULT 
OVERWRITES THE INPUT. 

MAXS ROW DIMENSION OF S 

IRS NUMBER OF ROWS IN s (irs.lf.maxs.and.irs.ge.2) 

JCS NUMBER OF COLUMNS IN S 

JSTART INDEX OF THE FIRST COLUMN TO BF. TRIANGULARIZED. JF 

JSTART. LT.l IT IS ASSUMED THAT JSTART=1, I.E. 

START TrIANGULARIZATION At COLUMN I. 

JSTOP INDEX OF LAST COLUMN TO Be TRIANGULARIZED. 

IF JSTOP. LT. JSTART. OR. JSTOP, GT. JCS THEN 
IF IRS. LE. JCS JSTOP IS SET EQUAL TO IRS-1 
IF IRS. GT. JCS JSTOP IS SET EQUAL TO JCS 
I.E. THE TRIANGULARIZATION IS COMPLETED AS FAR 
AS POSSIBLE 
V ( IRS) WORK VECTOR 

l 

COGNIZANT PERSONS: G. J .RlERMAM/M. W.NEAD (JPL. FEB. 1978) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION S(MAXS*JCs>* V(IRS> 

DOUBLE PRECISION SUM* DELTA 

ONE=l.oO 

ZERO=O.DO 

JSTT=JSTART 

JSTP=JSTOP 

IF (JSjT.LT. 1) JSTT=1 

IF (JSTP.GE.JSTT.AND. JSTP.LE.JCS) GO TO 5 
IF (IRS. LE. JCS) JSTP=IRS-1 
IF (IRS.GT.JCS) JSTP=JCS 

5 DO 40 J=JSTT*JSTP 
SUM=ZERO 
DO lo I=J*IRS 
v<D=s<I*J) 

S(l*J)=ZERO 

10 SUM=SUM+V(I)**2 

IF (SUM. LE. ZERO) GO TO 4q 

IF SUM=ZERO* COLUMN J IS ZERO AND THIS STEP OF THE 
ALGORITHM IS OMITTED 
SUM=DSQRT(SUM) 

IF (V<J> .GT.ZERO) sum=-sum 

S(J*J)=SUM 

V(J)=V(J)-SUM 


TDHHT010 
TDHHT020 
TDHHT030 
TDHHT040 
TDHHT050 
TDHHT060 
TDHHT070 
THHHT080 
TDHHT090 
TDHHTJ 00 
TDHHT110 
TDHHT120 
TDHHT130 
TDHHT140 
TDHHT150 
TDHHT160 
TDHHT170 
TDHHT1B0 
TDHHT190 
TDHHT200 
TDHHT210 
TDHHT220 
TDHHT230 
TDHHT240 
TDHHT250 
TDHHT260 
TDHHT270 
TDHHT2B0 
TDHHT290 
TDHHT300 
TDHHT310 
TDHHT320 
TDHHT330 
TPHHT340 
TDHHT350 
TDHHT360 
TDHHTV70 
TnHHT380 
TDHHT390 
TOHHT400 
TDHHTU10 
THHHT420 
TPHHT430 
TnHHT440 
TOHHT450 
TDHHT460 
THMHT470 
TOHHT480 
TDHHT490 
TDHHTSOO 
TPHHTS10 
TDHHT520 
TDHHT530 
TPHHT540 
TPHHTS50 
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SUM=0NE/CSUM*V( J) ) 

TDHHT560 


the householder transformation ts t=i-sum*v*v**t 

TPHHT570 


JP1=J+1 

TDHHTS80 


IF CjPl.GT.JCS> GO TO 40 

TDHHT590 


DO 30 K=JP1»JCS 

TDHHT600 


delta=zero 

TDHHT610 


DO 20 I=J*IRS 

TDHHT620 

20 

delta=delta+s ( 1 1 K > *V ( I ) 

TDHHT630 


DELTA=DELTA*SUM 

TDHHT640 


DO 30 I = J t IRS 

TDHHT650 

30 

S(I*K)=S(I»K)+DELTA*V(I> 

TDHHT660 

40 

continue 

T0HHT670 

TDHHT680 


RETURN 

TDHHT690 


END 

TDHHT700 
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page 1S 




SUBROUTINE THH (R>N»A*IA*M» SOS » NSTPT ) 

THIS SUBROUTINE PERFORMS A TR I ANGUL AR T Z AT I ON OF A RECTANGULAR 
MATRIX INTO A SINGLY-SUBSCRIPTED ARRAY BY APPLICATION OF 
HOUSEHOLDER ORTHONORMAL TRANSFORMATIONS. 

R(N*(N+3)/2> VECTOR STORED SQUARE ROOT INFORMATION MATRIX 

(LAST N LOCATIONS MAY CONTAIN A RIGHT HAND SIRE) 
N DIMENSION OF R MATRIX 

A(M»N+1) MEASUREMENT MATRIX 

IA ROW DIMENSION OF A 

M 


THHOOO10 

THH00020 

THH00030 

THH00040 

THH00050 

THH00060 

THHOOO^O 

THH00080 

THH00090 

THH00100 

THHOOllO 


NSTRT 


NUMBER OF ROWS OF A THAT ARE TO RE COMBINED WITH R THH00120 
(M.LE.IA) THH00130 
ACCUMULATED ROOT SUM OF SQUARES OF THE RESIDUALS THH00140 
SORT (Z-A*X (EST) **2) » INCLUDES A PRIORI THH00150 
SOS MUST BE INPUT AS A VARIARLE* NOT AS A THH00160 
NUMERICAL VALUE. IF INPUT S05.LT. Or MO SOS THH00170 
COMPUTATION OCCURS. THH00180 
FIRST COL OF THE INPUT A MATRIX THAT HAS A NONZERO THH00190 
ENTRY, IF NSTRT. Le.1» IT IS SET TO 1. THIS OPTION THH00200 


ENTRY, IF NSTRT. Le.1» IT IS SET TO 1. THIS OPTION THH00200 
IS CONVENIENT WHEN PACKING A PRIORI BY BATCHES AND THH00210 
THE A MATRIX HAS LEADING COLUMNS OF ZEROS. THH00220 

THH00230 

THH0024R 

ON ENTRY R CONTAINS A PRIORI SQUARE ROOT INFORMATION FILTER (SRIF)THHO0?5O 
ARRaY» AND ON EXIT IT CONTAINS THE A POSTERIORI (PACKED) ARRAY .THH00260 
ON ENTRY A CONTAINS OBSERVATIONS WHICH ARE DESTROYED BY THE THH00270 

INTERNAL COMPUTATIONS. THH00280 

ON ENTRY IF SOS IS .LT. ZERO * PROGRAM will ASSUME THERE IS NO THH00290 

RIGHT HAND SIDE DATA AND WILL NOT ALTER SOS OR USE LAST N THH00300 

LOCATIONS OF VECTOR R. THH00310 


EPS=-l,D-200 

ZERO=O.DO 

ONE=1.DO 

NSTART=NSTRT 


0 MACHINE DEPENDENT ACCURACY TERM 


IF (NSTART.LE.O) nstart=i 
NP1=N+1 

IF ( SOS, LT. ZERO) NP1=N 
KK=NSTART* (NSTART-1 ) /2 
DO 100 j=nstart*n 

KK=KK+j 
SUM=ZERO 
DO 20 1=1 fM 
20 SUM=SUM+A(I»J)**2 

IF (SUM, LE. ZERO) GO TO 100 

SUM=SUM+R(KK)**2 

SUM=DSQRT(SUM) 


ATIONS WHICH ARE DESTROYED BY THE THH00270 

THH00280 

*0 * PROGRAM wxlL ASSUME THERE IS NO THH00290 

D WILL NOT ALTER SOS OR USE LAST N THH00300 

THH003I 0 
THH00320 

cognizant persons g.j.bierman/n.hamata (Jpl» march 197 a) THH 00330 

THH00340 

IMPLICIT DOUBLE PRECISION (A-HrO-7) THH00350 

DIMENSION A(IA»1)»R(1) THH00360 

DOUBLE PRECISION SUM» ONE* BETA* DELTA THH00370 

THH00380 

;HINE DEPENDENT ACCURACY TERM THH00390 

THH00400 
THH00410 
THH00420 
THH00430 
THH00440 

13 NO. COLUMNS OF R THH00450 

0 NO COLS* = N IF SOS.LT.O THH00460 

THH00470 

0 J“TH STEP OF HOUSEHOLDER REDUCTION THH00480 

THH00490 
THH00500 
THH00510 
THH00>520 

fl IF J-TH COL. OF AiEQ.O GO TO STEP J+1THH00530 

THH00540 
THH00S50 


13 NO. COLUMNS OF R 

13 NO COLS* = N IF SOS.LT.O 

0 J-TH STEP OF HOUSEHOLDER REDUCTION 
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IF(R(KK) .GT.ZERO) SIJM=-SUM 

DELTA=R(KK>-SUM 

R(KK)=SUM 

JP1=J+1 

IF (JP1.GT.NP1) GO TO 105 

beta=sum*oelta 

IF (BETA. GT. EPS) GO TO 100 

beta=one/beta 

JJ=KK 

L=J 

C ** READY TO APPLY J-TH HOUSEHOLDER TRANS. 

DO 40 k=upi»npi 
JJ=JJ+L 
L=L+1 

SUM=DElTA*R(JJ) 

DO 30 I=1»M 

30 SUM=SUM+A(I» J)*A.(I»K> 

IF(SUM.EQ.ZERO) GO TO 40 
SUM=SUM*BETA 

BETA DIVIDE USED HERE TO AVOID OVERFLOW IN 
PROBLEMS WITH NEAR COLUMN COLLINEAPITY. IN THAT CASE 
COMMENT OUT LINE 630 AND CHANGE * TO / IN LINE 740 
R ( J J ) =R ( J J ) +5UM*DELTA 
DO 35 I=1»M 

35 A(I»K)=A(lFK)+SUM*A(If J) 

40 CONTINUE 
100 CONTINUE 

105 IF(SOS.LT.ZERO) RETURN 

CALCULATE SOS 

SUM=ZErO 
DO 110 1=1 f M 
110 SUM=SUm + A ( I * NP1 ) **2 
S0S=DSqRT ( SOS**2+SUM ) 

C 

RETURN 

END 


THH00560 

THH00570 

THH00580 

THH00590 

THH00600 

THH00610 

THH00620 

THH00630 

THH00640 

THH00650 

THH00660 

THH00670 

THH00680 

THH00690 

THH00700 

THH00710 

THH00720 

THH00730 

THH00740 

THH00750 

THH00760 

THH00770 

THH00780 

THH00790' 

THH00800 

THH00810 

THH00820 

THH00830 

THH00840 

THH00850 

THH00860 

THH00870 

THH00880 

THH00890 

THH00900 

THH00910 

THH00920 

THH00930 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


SUBROUTINE TTHH(R#rA#N) 

THIS SUBROUTINE COMBINES TWO SINGLE SUBSCRIPTED SRIF ARRAYS 
USING HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS 


R(N*(N+l)/2) 

RA(N*tN+l)/2) 

N 


INPUT VECTOR STORED UPPER TRIANGULAR MATRIX# 
RESULT IS IN R 

THE SECOND INPUT VECTOR STORED UPPER TRIANGULAR 
MATRIX. THIS MATRIX IS DESTROYED BY THE 
COMPUTATION 

DIMENSION OF THE ESTIMATED PARAMETER VECTOR. 

A NEGATIVE VALUE FOR N IS USED TO NOTE THAT 
R AND RA HAVE RT. HAND SIDES INCLUDED AND 
HAVE DIM=ARS(N)*(ABS(N)+3)/2i 

ON EXIT RA IS CHANGED AND R CONTAINS THE RESULTING SRlF ARRAY 

COGNIZANT PERSONS G. J.BlERMAN/M.W.NEAD (JPL» JAN. 1976) 

IMPLICIT DOUBLE PRECISION! A-H.O-Z) 

DIMENSION RA(1>» R{1) 

DOUBLE PRECISION SUM 0 FOR USE IN SINGLE PRECISION VERSION 

ZERO=0. 

one=i. 

NP1-N 

IF (N.eT.O) GO TO 10 
N=-N 


0 I J(START) 

Q J-TH STEP OF HOUSEHOLDER REDUCTION 


NPl=N+l 
10 IJS=1 
KK=0 

DO 100 J=1»N 
KK=KK+J 
SUM=R<KK>**2 
DO 2o I=IJS»KK 
20 SUM=SuM+RA(I>**2 

IF (SUM. LE. ZERO) GO TO 100 
SUM=SqRT(SUM) 

IF (R(KK) .GT.ZERO) SUM=-SUM 

DELTA=R(KK)-SUM 

R(KK)=SUM 

BETA=0NE/ ( SUM*DELTA ) 

JJ=KK 

L=J 

JP1=J+1 

IKS=KK+1 

C * * * J-TH HOUSEHOLDER TRANS. DEFINED 

C 40 LOOP APPLIES TRANSFORM* TO COLS. J+l TO NPl 

DO 40 K=JP1#NP1 
JJ=JJ+L 
L=L+1 
IK=IKS 

SUM=DELTA*R(JJ) 

DO 30 I=IJS»KK 
SUM=SuM+RA ( IK) *RA ( I ) 


TTHH0010 

TTHH0020 

TTHH0030 

TTHH0040 

TTHH0050 

TTHH0060 

TTHH0070 

TTHH0080 

TTHH0090 

TTHH0100 

TTHH0110 

TTHH0120 

TTHH0130 

TTHH0140 

TTHH0150 

TTHH0160 

TTHH0170 

TTHH0180 

TTHH0190 

TTHH0200 

TTHH0210 

TTHH0220 

TTHH0230 

TTHH0240 

TTHH0250 

TTHH0260 

TTHH0270 

TTHH0280 

TTHH0290 

TTHH03OO 

TTHH0310 

TTHH0320 

TTHH033O 

TTHH0340 

TTHH0350 

TTHH0360 

TTHH0370 

TTHH0380 

TTHH0390 

TTHH0400 

TTHH0410 

TTHH0420 

TTHH0430 

TTHH0440 

TTHH0450 

TTHH0460 

TTHH0470 

TTHH0480 

TTHH0490 

TTHH05OO 

TTHH051O 

TTHHO820 

TTHH0530 

TTHH0540 

TTHH0550 
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30 

IK=1K+1 

TTHH0560 


IF (SuM'.EQ.ZERO) go TO 40 

TTHH0^70 


SUM=SuM*BETA 

TTHH0«5Q0 


R ( J J ) =R ( J J ) +SUM+DEL.TA 

TTHH0590 


IK-IKS 

TTHH0600 


DO 35 I=IJS»KK 

TTHH0610 


RA(IK)=RA(IK)+SUM*RA(I) 

TTHH0620 

35 

IK=IK+1 

TTHH0630 

40 

IKS=I«S+K 

TTHH0640 

100 

IJS=KK+1 

TTHH0650 

TTHH0660 


RETURN 

TTHH0670 


END 

TTHH0680 
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ORIGINAL PAGE IS 
OF POOR QUALITY! 


SUBROUTINE TWOMAT (A»N»LENrCAR»TEXT»NCHAR»NAMES) 

to display a vector stored upper triangular matrix in a 

TWO-DIMENSIONAL TRIANGULAR FORMAT 

A (N*(N+l)/2) VECTOR CONTAINING UPPER TRIANGULAR MATRIX (Dp) 

N DIMENSION OF Matrix (I) 

LEN NUMBER OF COLUMNS TO BE PRINTED* 7 OR 12 (I) 

CAR (N) PARAMETER NAMES (I) 

TEXT( ) An ARRAY of FTELHATA CHARACTERS TO RE PRINTED AS 

A TITLE PRECEDING TH^ MATRIX 

NCHAR NUMrER OF CHARACTERS' INCLUDING SPACES* THAT 

ARE to BE PRINTED IN TEXT ( ) 

ABS(NCHAR) .LE.114. NCHAP NEGATIVE IS USED 
TO AVOID SKIPPING TO A NEW PAGE TO START 
PRINTING 

names true to print paramfTer names 

COGNIZANT PERSON! M.W.NEAD (JPL» OCT. 1977) 

PARAMETER J12=12» J7=7 
DOUBLE PRECISION A (N) 

INTEGER CARCN)* TEXtU)' L(Jl2)» LIST(Jl2) 

LOGICAL NAMES 

INTEGER V (4 ) *VFMT(jl2) *V7MT(J7) *V1 rMT<J12) 

DATA v/*(2X**»»A6rlX**#» * * *E10 .5) »/ * < V12MT ( I ) * 1=1 * 12) 

1 /»12* , ' 10 X* 11 * * ^OX'lO* * *30X»9' »»040X*e* * '050Xf7* , 

2 *060X»6* * *070X*5» * *080X*4* * »090X*3* * *100X*2» * f llOX.l*/* 

1 V7MT/«7» , *017X*6' * *034X*5» * * 051X * 4 * • * 068X *3» * '085X*2* * *102X,1*/ 
DATA KON7/*D17.85 f /> KON12/ * ElO . 5) * / 

M1*m2 ROW LIMITS FOR EACH PRINT SEQUENCE ’ 

N1»M2 COL LIMITS FOR EACH LINE OF PRINT 

L( I ) LOC OF EACH COLUMN IN A ROW 

KT ROW COUNTER 

* * * * * INITIALIZE COUNTERS 

IF (LEN.EQ.J0) GO TO 5 
IF (LEN.EQ.7) GO TO 1 
IF (LEN.EQ.12) GO TO 2 
WRITE (6*230) LEN 
LEN=12 
GO TO 2 

1 V(4)=K0N7! JO=75 J0Ml=J0-ll JOP1=JO+1 » 

1 REPEAT 1=1* JO) VFMT(I)=V7MT(I) 

GO TO 5 

2 V ( 4 ) =K0N12 ) J0=l2) J0M1=J0-1) JOPl=jO+l) 

1 REPEAT I=1*J0) VFMT(I)=V12MT(I) 

5 Ml=l 
M2=J0 
Nl=l 
KT=0 

V (2 ) = ' A6* IX* * 

IF (.NOT. NAMES) V(2)='I5*2X* 


TWOMOOIO 

TWOM0020 

TWOM0030 

TWOM0O4O 

TWOM0050 

TWOM0060 

TWOM0070 

TWOM0080 

TWOM0090 

TWOMOIOO 

TWOM0110 

TWOM0120 

TWOM013O 

TWOM0140 

TWOM0150 

TWOM0160 

TWOM0170 

TWOM0180 

TWOM0190 

TWOM0200 

TWOM0210 

TWOM0220 

TWOM0230 

TWOM0240 

TWOM0250 

TWOM0260 

TWOM0270 

TWOM0280 

TWOM0290 

TWOM0300 

TWOM031O 

TWOM0320 

TWOM0330 

TWOM0340 

TWOM0350 

TWOM0360 

TWOM0370 

TWOM0380 

TWOM0390 

TWOM0400 

TWOM0410 

TWOM0420 

TWOM0430 

TWOM0440 

TWOM0450 

TWOM0460 

TWOM0470 

TWOM0480 

TWOM0490 

TWOM0500 

TWOM0510 

TWOM0520 

TWOM0530 

TWOM0540 

TWOM0550 
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NC=IABS<NCHAR)/6 

IF (MOD (NCHAR.6) »NE.O) NC=NC+1 

IF (NCHAR.GE.O) WRITE (6.200) (TEXT ( I ) » 1=1 »NC) 

IF (NCHAR.LT. 0) WRITE (6.205) (TEXT ( I ) » 1=1 »NC ) 

10 IF (M2.GT.N) M2=N 

IF (.NOT. NAMES) GO TO 20 

IF <LEN*EQ.7) WRITE (6.210) (CAR ( I ) . I=N1 .M2) 

IF (LEN'EQ.12) WRITE (6»2il) (CAR ( I ) ' I=Nl »M2) 

GO TO 40 
20 M=Nl 

L2=M2-N1+1 
DO 30 1=1 »L2 
LIST ( I )=M 
30 M=M+1 

IF (LEN«EQ.7) WRITE (6.220) (LIST ( I ) » 1=1 .L2) 

IF (LEN.EQ.12) WRITE (6»22l) (LIST ( I ) . 1=1 »L2) 

40 CONTINUE 
C ****** 

DO 190 IC=M1»M2 
K=1 

IF (iC.LE. (KT+JO) ) GO TO 60 
JJ=0 

DO 5o J=1»IC 
JJ=JJ+J 
L(K)=JJ 
I1=IC"KT*J0 

IF (ll.EQ.JO) 60 TO 90 
GO TO 70 
CONTINUE 

11=1 

LOO zL 00+1 
CONTINUE 
DO 8o I=Il»JOMl 
K=K + 1 

II=I+KT*JO 

LOO=L(K-1)+II Q OBTAIN COL INDEX FOR ROW 

CONTINUE 

I2=MiNO(JOPl» (M2+l-KT*J0) )-Il 
V ( 3) zVFMT (ID 
IF (.NOT. NAMES) GO TO 180 
WRITE (6. V) CAR(IC)»(A(L(I))0=1.I2) 

GO To 190 

WRITE (6.V) IC » ( A (L ( I ) ) . 1=1 . 12) 

CONTINUE 

IF (M2.EQ.N) RETURN 
N1=M2+1 
M2ZM2+J0 
KT=KT+1 

IF (NCHAR.GE.O) WRITE (6.201) (TEXT( I ) . 1=1 .NC) 

IF (NCHAR.LT.O) WRITE (6.206) (TEXT ( I ) . 1=1 #NC ) 

GO TO 10 

200 FORMAT (1H1.2X.21A6) Q TITLE 

205 FORMAT (1H0.2X.21A6) 0 TITLE 


50 


60 

C 


70 


80 

'90 


180 

190 


TWOM0560 

TWOM0570 

TWOM0580 

TWOM0590 

TWOM0600 

TWOM0610 

TWOM0620 

TWOM0630 

TWOM0640 

TWOM0650 

TWOM0660 

TWOM0670 

TW0M0680 

TWOM0690 

TWOM0700 

TWOM0710 

TWOM0720 

TWOM0730 

TW0M0740 

TWOM0750 

TWOM0760 

TW0M0770 

TWOM0780 

TW0M0790 

TWOMOSOO 

TWOM0810 
TWOM0820 
TWOM0830 
TWOM0840 
TWOM0850 
TWOM0860 
TWOM0R70 
TW0M0880 
TWOM0890 
TW0M0900 
TW0M0910 
TW0M0920 
TW0M0930 
TWOM0940 
TWOM0950 
TWOM0960 
TWOM0970 
TWOM0980 
TWOM0990 
TWOMIOOO 
TWOMIOIO 
TWOM1020 
TWOM1030 
TWOM1040 
TWOM1050 
TWOM1060 
TWOM1070 
TWOM1080 
TWOM1090 
TWOMllOO 
TWOM1110 
TWOMl 120 


124 



™>traWAL PAGBlS 


201 

FORMAT 

( lt+1 * 2X t » (CONTINUE) 

* »19A6) 

R TITLE 

TW0M1130 

206 

FORMAT 

(1H0*2Xf * (CONTINUE) 

* » 19A6) 

Q TITLE 

TW0M1140 

210 

FORMAT 

(1H0#5Xf7(HX»A6)) 

Q HORIZONTAL 

names 

TWOM1150 

220 

FORMAT 

(1H0»3X>7(11X»I6)) 



TW0M1160 

211 

FORMAT 

(1H0f5Xf12(4x>A6)) 

ra horizontal 

NAMES 

TWOM1170 

221 

FORMAT 

(1H0»3X»12(4x»I6)> 



TW0M1180 

230 

FORMAT 

(1H0»20Xf*TWOMAT CALLED WITH LENGTH = 

* » 13) t 

TWOM1190 

C 

END 

- 



TWOM1200 

TWOM1210 
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SUBROUTINE TZERO (RfNrlSrlF) - TZEROOOO 

TZEROOlO 

TO ZERO OUT ROWS IS (ISTART) TO IF ( IFINAL) OF A VECTOR .TZER0020 

STORED UPPER TRIANGULAR MATRIX TZER0030 

TZERO040 

R(N*(N+l)/2) INPUT VECTOR STORED UPPER TRIANGULAR MATRIX TZER0050 

M DIMENSION OF R . TZER0060 

IS FIRST row of R THAT IS TO be Set to ZERO TZER0070 

IR LAST ROW OF R that Is TO BE SET TO ZERO TZERO080 

TZER0090 

COGNIZANT persons: G.J.bIERMAN/C«F. PETERS <JPL> NOV. 1975) TZEROIOO 

TZERO110 

IMPLICIT DOUBLE PRECISION (A-H»0-Z) TZER0120 


DIMENSION R ( 1 ) TZER0130 

C T7ER0140 

ZERO=O.DO TZERO150 

IJS=lS*(IS-l)/2 TZERO160 

DO lo I=IS»IF TZER0170 

IJS=lJS+I TZERO180 

IJ=IJS TZER0190 

DO 10 J=I*N TZER0200 

R(IJ)=ZERO TZERO210 

IJ=IJ+J TZER0220 

10 CONTINUE TZERO230 

C TZER0240 

RETURN TZER0250 

END TZER0260 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


SUBROUTINE UDCOL (U*N* KS* NCOLOR * V *EM* Q) 

COLORED NOISE UPDATING OP THF U-D covariance FACTORS* I.E. 
U*D* (U**T) -OUTPUT=PHI*U*D* (U**T)* (PHI**T) +Q 
PHI=DlAGtO(KS“l) >EM(1) ♦ . »EM(NCOLOR) .0(N-(KS-1+NC0L0 r)) ) 
Q=DIAG(0 (KS-1) .0(1) » . . . *Q(NCOLOR) .0(N-(kS-1+NC0LOR) ) ) 

O(K) Is A VECTOR OF ZEROS 

The ALGORITHM USED IS THE BIERMAN-ThORNTON ONE COMPONENT 
AT-A-TlME UPDATE • CF • BIERMAN ^FACTORIZATION METHOD 
FOR DISCRETE SEQUENTIAL ESTIMATIONS* ACADEMIC PRESS C1977) 
PP. 147-148 


U(N*(N+l)/2) 


KS 

NCOLOR 

V(KS-l+NCOLOR> 
EM (NCOLOR) 

Q( NCOLOR) 


INPUT U-D VECTOR STORED COVARIANCR FACTORS. 

THF COLORED NOISE UPDATE RESULT RESIDES 
IN U ON OUTPUT 

FILTER DIMENSION. IF THE LAST COLUMN OF U 

HOUSES THF FILTER ESTIMATFS* THEN 

N=NUMBER FILTER VARIABLES + 1 

THE LOCATION OF THE FIRST COLORED NOISE TERM 

<KS.GE,1.AND.KS.LE»N) 

The NUMBER OF COLORED NOISE TERMS (NCOLOR .6E. 1 > 

work vector 

INPUT VECTOR OF fcoLORED NOISE MAPPING TERMS 

(UNALTERED BY PROGRAM) 

input vector of process noise variances 

(UNALTERED BY PROGRAM) 


SUBROUTINE REQUIRED: RANK! 

COGNIZANT PERSON! G.J. BIERMAN (JPL* JAN. 1978) 
DOUBLE PRECISION TMP*S 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION U(l) »V(1) ,EM(1> *Q(1) 

****** INITIALIZATION 
NM1=N-1 
KSM1=KS-1 
JJ0LD=KS*KSMl/2 
KOL=KSMl 
****** 


00 50 K = 1»NC0L0R 
K0LM1=K0L 
K0L=K0L+1 
JJ=JJ0LD+K0L 
TMP=U(JJ)*EM(K) 
C=Q(«)*U(JJ) 
S=TMp*EM(K)+Q(K) 
U(JJ)=S 


QD(J) update 


IF (KOL.GE.N) go TO 20 
IJ=JJ 

DO lo J=K0L*NM1 
IJ=IJ+J 


uncoLniO 

UDCOL020 

UDC0L030 

UDCOLO40 

UDC0L050 

UncOLOGO 

UDCOL070 

UDC0L080 

UDCOL090 

UncOLlOO 

UDCOL110 

UDCOL120 

UDCOL130 

UDCOL140 

UDCOL150 

UnC0Ll60 

UDC0L170 

UDC0L180 

UDCOL190 

UDC0L200 

UDC0L210 

UDC0L220 

UDC0L230 

UDC0L240 

UDC0L250 

UDC0L260 

UDC0L270 

UDCOL280 

UDC0L290 

UDC0L300 

UDC0L310 

UDC0L320 

UDCOL330 

UDC0L340 

UDCOL350 

UDCOL360 

UDCOL370 

UDC0L380 

UDC0L390 

UDCOL400 

U0C0L410 

UDCOL420 

UDC0L430 

UDC0L440 

UOCOL450 

UDC0L460 

UDC0L470 

UDC0L480 

UDC0L490 

UDC0L500 

UDCOL510 

UDC0L520 

UDCOL530 

UDC0L540 

UDC0L550 



10 U<IJ)=U(IJ)*EM(K) 13 UPDATING ROW KOL ENTRIES 

C 

20 IF (JJ.EQ.1) GO TO 50 0 (WH^N KS=1» N-l) 

IF (S»LE.0.D0) GO TO 30 

TMP=TMP/S 0 TMP=EM<K)*D(KOL) -OLD/D (KOL) -NEW 

C=c/S Q C=QlK>*D(KOL)-OLD/D(KOLJ-NEW 

30 DO 40 I = 1 f KOLMl 

V(l)=U(JJOLD+I) 

40 U( JJ0LD+I)=TMP*V(I) 

IF (KOLMl. GT.l) GO TO 45 
U(1)=U<1)+C*V<1)**2 
GO TO 50 

45 CALL RANK1(U*U>K0LM1»C,V) 

50 JJ0LD=JJ 

C 

RETURN 

END 


UDC0L560 

UDC0L570 

UDCOL580 

UDCOL590 

UDC0L600 

UDCOL610 

UDCOL620 

UDCOL630 

UDCOL640 

UDCOL650 

UDCOL660 

UDC0L670 

UDCOL680 

■UDCOLf.90 

UDC0L700 

UDCOL710- 

UDCOL720 
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SUBROUTINE UDMEAS (U*N*R*AfF*G* ALPHA ) 

COMPUTES ESTIMATE AMD U"D MEASUREMENT UPDATED 
COVARIANCE* P=U*D*U**T 

*** INPUTS *** 

U UPPER TRIANGULAR MATRIX t WJTh p ELEMENTS STORED AS THE 

DIAGONAL. U IS VECTOR STORED AND CORRESPONDS TO THE 
A PRIORI COVARIANCE* IF STATE ESTIMATES ARE COMPUTED* 
THE LAST COLUMN OF U CONTAINS X. 

N DIMENSION OF THE STATE ESTIMATE. N.GT.'l 

R MEASUREMENT VARIANCE 

A VECTOR OF MEASUREMENT COEFFICIENTS* IF DATA THEN ACN+D 

ALPHA IF ALPHA LESS THAN ZERO NO ESTIMATES ARE COMPUTED 
(AND X AND 2 NEED NOT BE INCLUDED) 

*** OUTPUTS *** 

U UPDATED* VECTOR STORED FACTORS AND ESTIMATE AND 

U( (N+l) (N+2)/2) CONTAINS (Z-a**T*X) 

ALPHA INNOVATIONS VARIANCE OF THE MEASUREMENT RESIDUAL 
G VECTOR OF UNWEIGHTED KALMAN GAINS. THE KALMAN 

GAIN K IS FQUAL TO G/ALPHA 
F CONTAINS IJ**T*A AND (Z-A**T*X) /ALPHA 

ONE CAN HAVE F OvrRWRlTE A TO SAVE STORAGE 

cognizant persons: g.j. bierman/m.w. nfad cjpl* feb.1978) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION U(l)» A ( 1 ) . F(l), G(l) 

DOUBLE PRECISION SUM* BETA* GAMMA 
LOGICAL IEST 

ZEROrO.DO 
IEST=. FALSE. 

one=i.do 

NP1=N+1 

NP2=N+2 

NT0T=N*NPl/2 

IF (ALPHA.LT, ZERO) GO TO 3 
SUM=A(NP1) 

DO 1 J=1*N 

SUM=SUM-A ( J ) *U ( NTOT+J ) 

U(NT0T+NP1)=SUM 0 Z=Z-A**T*X 

TEST=.TRUE. 


o Z=Z-A**T*X 


3 JJN=NT0T 
DO 10 L=2*N 
J=NP2-L 

SUm=A(J) 

JMl-J-1 

DO 5 K=1»JM1 


UDMEA010 

UDMEA020 

UDMEA030 

UDMEA040 

UDMEA050 

UDMEA060 

UDMEA070 

UDMEA080 

UDMEA090 

UDMEA100 

UDMEA110 

UDMEA120 

UDMEA130 

=ZUDMEA140 

UDMEA150 

UDMEA160 

UDMEA170 

UDMEA180 

UDMEA190 

UDMEA200 

UDMEA210 

UDMEA220 

UDMEA230 

UDMEA240 

UDMEA250 

UDMEA260 

UDMEA270 

UDMEA280 

UDMEA290 

UDMEA300 

UDMEA310 

UDMEA320 

UDMEA330 

UDMEA340 

UDMEA350 

IJDMEA360 

UDMEA370 

UDMEA380 

UDMEA390 

UDMEA400 

UDMEA410 

UDMEA420 

UDMEA430 

UDMEA440 

UDMEA450 

UDMEA460 

UDMEA470 

UDMEA480 

UDMEA490 

UDMEA500 

UDMEA510 

UDMEA520 

UDMEA530 

UDMEA540 

UDMEA550 



o o o 


5 SUM=SUN|+U( JJ+)0*A<K) 

F(j)-=SUM 
G<J)=SuM*U( JJN) 

10 JJN=JJ 
F ( 1 ) = A (1 ) 

_G { 1 ) =U ( 1 ) *F ( 1 ) 

C - F=U**T*A AND G=D*(U**T*A) 

c 

sum=r+g(1)*f(d q sum c i ) 

GAMMA=0 0 FOR R=0 CASE 

IF <5Um*GT.ZER0) GAMMA=ONE/SUM 0 FOR R=0 CASE 

IF (F(l) .NE.ZERO) U ( 1 )=U( 1 ) *R*GAMMA 0 D(l) 

C • . 

KJ=2 

DO 20 J=2»N 
RETA=SUM 
TEMP=G(J) 

SUm=SUM+TEMP*F(J) 
p=_F(J>*GAMMA 
JMl=J-l 
DO 15 K-l t JMl 
S=U(KJ) 

U(KJ)=S+P*G(K) 

G(K)=G(K)+TEMP*S 
15 ' - KJ=KJ+1 

IF ( TEMP. EQ. ZERO) GO TO 20 
• GAmMA=ONE/SUM 

u<kJ)=u(kj)*beta*gamma 

20 KJ=KJ+i 

, alpha=sum 

* EQN. NOS. REPER TO BIERMAN’S 1975 CDC PAPER* PP. 337-346, 

IF (.NOT.IEST) RETURN 
F(NPl)=U{NTOT+NPl)*GAMMA 
DO 30 J=1»N 

30 U CNTOT+J) =U (NTOT+J) +G ( J) *F CNP1) 

RETURN 
END 


0 BETA=SUM<J-1) 

Q SUM(J) 

0 P=-F(J)*(1/SUM(J-1)> EQM(21) 


0 EQN(22) 

0 Eon (23) 

' B FOR R=0 CASE 
Q GAMMA=1/SUM(J) 

0 D ( J)‘ EON ( 19) 


UDMEA560 

UDMEAS70 

UDMEAS80 

UDMEA590 

UDMEA600 

UDMEA610 

UDMEA620 

UDMEA630 

UDMEA640 

UDMEA650 

UDMEA660 

UDMEA670 

UDMEA680 

UDMEA690 

UDMEA700 

UDMEA710 

UDMEA720 

UDMEA730 

UDMEA740 

UDMEA750 

UDMEA760 ’ 

UDMEA770 

UDMEA780 

UDMEA790 

UDMEA800 

UnMEABlO 

UDMEA820 

UDMEA830 

UDMEA840 

UDMEA850 

UDMEA860 

UDMEA870 

(JDMEA880 

UDMEA890 

UDMEA900 

UHMEA910 

UDMEAP20 

UDMEA930 

UDMEA940 

UDMEA950 
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c 

c 


c 


SUBROUTINE UD2C0V (ijIN*P0UT*N) 

to obtain a covariance from its u-d factorization, both matrices 
are vector stored and the output covariance can overwrite the 

INPUT U“D array. UIN=U-D IS RELATED To POUT VIA POUT=UDU(**T) 

UIN(N*{N+l)/2) INPUT u-d factors* vector stored with the d 
Entries stored on the diagonal of uin 
P 0UT<N*(N+l>/2) output COVARIANCE* VECTOR STORED. 

(P0UT=UIN IS PERMITTED) 

N DIMENSION OF THE MATRICES INVOLVED* N.GT.l 

COGNIZANT PERSONS: G.U.BIERMAN/M.W.nEAD IJPL* FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

.DIMENSION UIN(l) * POUT(D 

POUTfl)=UIN(l) 

JJ=1 

DO 20 J=2*N 

UJL— jJ 0 

JJ-Jj+J 

POUT(JJ)=UIN(JJ) 

S=P0UT(JJ) 

11=0 
JMl=j-l 
DO 2o I=1*JM1 
II=II+I 

ALpHA=S*UIN(JJL+I) P jJL+I=(I*J) 

IKsII 

DO 10 K=I*JM1 

P0UT(IK)=P0UTCIK)+ALPHA*UIN(JJL+K) ra JJL+K-IK * J) 

10 IK=IK+K 

20 POUT(JJL+I)=ALPHA 

RETURN 

END 


UD2C0010 

UD2C0020 

UD2C0030 

UD2C0040 

UD2C0050 

UD2C0060 

UD2C0070 

UD2C0080 

Un2C0090 

UD2C0100 

UD2C0110 

UD2C0120 

UD2C0130 

UD2C0140 

UD2C0150 

UD2C0160 

UD2C0170 

UD2C0180 

UD2C0190 

UD2C0200 

UD2C0210 

UD2C0220 

UD2C0230 

U02C0240 

UD2C0250 

UD2C0260 

UD2C0270 

UD2C0280 

UD2C0290 

UP2CO30O 

UD2CO310 

UD2C0320 

UD2C0330 

UD2CO340 

UD2CO350 

UP2C0360 

UD2C0370 

UD2C0380 
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SUBROUTINE UD2SIG(U*N*SIG*TEXT*NCT) 


COMPUTE STANDARD DEVIATIONS (SIGMAS) FROM U-D COVARIANCE FACTORS 


U(N*(N+l)/2) 


N 

SIG(N) 
TEXT( ) 

NCT 


INPUT VECTOR STORED ARRAY CONTAINING THE U-D 

FACTORS. THF D (DIAGONAL) ELEMENTS ARE STORED 

ON THE DIAGONAL 

U MATRIX DIMENSION* N.GT.l 

VECTOR OF OUTPUT STANDARD DEVIATIONS 

ARRAY OF FIELDATA CHARACTERS TO RE PRINTED 

PRECEDING THE VECTOR OF SIGMAS 

NUMBER OF CHARACTERS IN TFXT* 0oLE.NCT.LE.126 

IF NCT=0 * NO SIGMAS ARE PRINTED 


COGNIZANT PERSONS! G. J.BIERMAN/M.W.NEAD (JPL* FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

INTEGER TEXT Cl) 

DIMENSION U( 1 > » SIG(1) 

JJ=1 

SIG ( 1 ) =U ( 1 ) 

, DO 10 J=2*N 

, JJL-JJ S) (J-1.J-1) 

JJsJj+J 
. S=U(JJ) 

SIG(J)=S 
,JM1=J-1 
* DO 10 1=1* JM1 

10 . SIg(I)=SIG(I)+S*U(JJL+I)**2 
WE NOW HAVE VARIANCES 
DO 20 J=1*N 

20 , SIG(J)=SQRT(SIG(J) ) 

IF (NCT* EG. 0) GO TO 30 
. NC-NCT/6 

IF ( MOD ( NC * 6 ) « NE * 0 ) NC=NC+1 
WRITE (6*40) (TEXT ( I ) * 1=1 * NC) 

WRITE (6*50) (SIG(I) *I=1,N) 

30 RETURN 

40 FORMAT ( 1H0 *2X *2lA6) 

50 FORMAT < 1H0 * (6018.10 ) ) 

END 


UD2SI010 

UD2SI020 

UD2SI030 

UD2SI040 

UD2SI050 

UD2SI060 

UD2SI070 

UD2SI080 

UD2SI090 

UD2SI100 

UD2SI110 

UD2SI120 

UD2SI130 

UD2SI140 

UD2SI150 

UD2SI160 

UD2SI170 

UD2SI1S0 

UD2SI190 

UD2SI200 

U02SI210 

UD2SI220 

UD2SI230 

UD2SI240 

UD2SI250 

UD2SI260 

UD2SI270 

UD2SI280 

UD2SI290 

UD2SI300 

UD2SI310 

UD2SI320 

UD2SI330 

UD2SI340 

UD2SI350 

UD2SI360 

UD2SI370 

UD2SI380 

UD2SI390 

UD2SI400 

UD2SI410 

UD2SI420 

UD2SI430 

UD2SI440 

UD2SI450 
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C 

C 


C 

C 


C 


C 

C 


SUBROUTINE UTlNV(RIN»N.ROUT) 

TO INVERT AN UPPER TRIANGULAR VECTOR STORED MATRIX AND STORE 
THE RESULT IN VECTOR FORM. THE ALGORITHM IS SO ARRANGED THAT 
THE RESULT CAN OVERWRITE THE INPUT. 

IN ADDITION TO SOLVE RX=Z. SET RlN<M*<N+l)/2+l)=Ztl) . ETC . . 
AND SET RIN( (N+l ) * (N+2) /2) =-l . CALL THE SUBROUTINE USING N+l 
INSTEAD OF N. ON RETURN THE FIRST N FNTRTES OF COLUMN N+l 
WILL CONTAIN X, 

RIN(N*(N+l)/2> INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
N MATRIX DIMENSION 

ROUT < N*< N+l >/2} OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
INVERSE 

COGNIZANT PERSONS! G. J.BIERMAN/M.W.NEAD ( JPL. JAN. 1978) 

DOUBLE PRECISION RlN(l). ROUT(l), ZERO. DINV. ONE. SUM 

ZERO=O.DO 

ONE=1.DO 

IF (RlNd> .NE.ZERO) GO TO 5 
J=1 

WRITE (6.100) J.J 
RETURN 

5 ROUT(l)=ONE/RIN(l) 

JJ=1 

DO 20 J=2.N 
JJOLD=JJ 
J J= J J+ J 

IF (RIN(JJ) .NE.ZERO) GO TO 10 

WRITE (6.100) J.J 

RETURN 

10 DINV=ONE/RIN(JJ) 

ROUT { JJ)=DINV 
11 = 0 
IK=1 
JM1=J-1 
DO 2o 1=1. JM1 
II=II+I 
IKrII 

sum=zero 

DO 15 K=I . JMl 

SUM=SUM+ROUT ( IK ) *RIN ( JJOLD+K ) 

15 IK=IK+K 

20 ROuT(JJ0LD+l)=-SUM*DTNV 

RETURN 


UTINV010 

UTINV020 

UTINV030 

UTINV040 

UTINV050 

UTINV060 

UTINV070 

UTINV080 

UTINV090 

UTINV100 

UTINV110 

UTINV120 

UTINV130 

UTINV140 

UTIMV150 

UTINV160 

UTINV170 

UTINV180 

UTINVI90 

UTINV200 

UTINV210 

UTINV220 

UTINV230 

UTINV240 

UTINV250 

UTINV260 

UTINV270 

UTIMV280 

UTINV290 

UTINV300 

UTIMV310 

UTINV320 

UTINV330 

UTJNV340 

UTINV350 

UTINV360 

IJTINV370 

UTINV380 

UTINV390 

UTINV400 

UTINV410 

UTINV420 

UTIMV430 

UTINV440 

UTINV450 

UTINV460 

UTINV470 

UTINV480 

UTINV490 

UTINV500 

UTINV510 

UTINVS20 


UTINV530 

loo format (iho.iox.'* * * matrix inverse computed only up to but mot utinv540 
IINCLUDiNG COLUMN’ .14.' * * * MATRIX DIAGONAL '»I4.» IS ZERO * * **UTINV550 
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UTINV560 

UTINV570 

UTINV580 
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C 


SUBROUTINE UTIROW (RINrNrROUT fNRY) 

TO COMPUTE THE INVERSE OF AN UPPER TRIANGULAR (VECTOR STORED) 
MATRIX WHEN THE LOWER PORTION OF THF INVERSE IS GlVFN 

ON input: 

RX RXY * * RX RXY 

Rin= rout= where r= 

* * o py**-i o ry 


on output: rin is unchanged and Rout=r**-i 

the RESULT can OVER-WRITE THE INPUT (I.F. RTNrROUT) 


RlN<N*(N+l)/2) 

N 

R0UT(N*(N+l)/2) 


NRY 


INPUT VECTOR STORED TRIANGULA P MATRIX 
THE BOTTOM NRY POV/S ARE IGNORED 

matrix dimension 

OUTPUT VECTOR STORED MATRIX. ON INPUT THE 

bottom nry rows contain the lower portion 
of R**-l. on output ROUT=R**-l 
DIMENSION of LOWER (ALREADY INVERTED) 

triangular r. if nRy=o» ordinary matrix 
INVERSION RESULTS. 


COGNIZANT PERSONS: G.U.BIERMAN/M.W*NEAD (JPL MARCH 1977) 

double precision rin(1)» rout(d* sum. zerOf onef dinv 
data one/i. do/» ZERO/0. do/ 


INITIALIZATION 


c 


c 


NR=N*(N+l)/2 
istrt=n-nry 
irlst=istrt+i 
II=ISTRT* lRLST/2 
DO 40 lROW=ISTRT»lr-l 
IF (RlN(II) .NE*ZERO) go 
write (6f50) IROW 
return 

10 DINV=ONE/RINlin 

ROUT ( I I ) = dinv 

KJS=NR+IROW 

IKS=lI+IROW 


0 No. ELEMENTS In r 
ra FIRST row to be inverted 

0 IRLST=PREVI0US IROW INDEX 
0 II=DlAGONAL 

TO 10 


0 KU(START) 
0 IK (START) 


IF (iRLST.GT.N) GO TO 35 
DO 30 J=N»IRLSTf-1 
kjs=kjs-j 
sum=zero 

IK=IKS 

KJ=KJS 


DO 20 K=IRLSTfJ 

kj=kj+i 

SUM=SUM+RIN ( IK) *ROUT (KJ) 


UTIPOOOO 
UTIRCOlO 
UTIR0020 
SJTIR0030 
IJTIR0040 
UTIR0050 
UTIR0060 
UTIR0070 
UTIROO80 
IJTIRO090 
UTIRO) 00 
UTIRO110 
UTIR0120 
UTIRO130 
UTIR0140 
UTIRO) 50 
UTIR0160 
UTIRO) 70 
UTTR0180 
UTIR0190 
UTIR0200 
IJTIROPlO 
UTIR0220 
UTIR0230 
UTIRO240 
UTIR0250 
UTIR0260 
UTIR0270 
UTIR0280 
UTIRO290 
UTIP0300 
UTIR0310 
UTIRO320 
UTIRO330 
UTIP0340 
UTIR0350 
UTIR0360 
UTIR0370 
UTIR0380 
UTIR0390 
UTIR0400 
UT1R0410 
UTIR0420 
UTIRO430 
UTIR0440 
UTIRO450 
UTIR0460 
UTIRO470 
UTIR0480 
UTIR0490 
UTIRO50O 
UTIRO510 
UTIR0520 
UTIR0530 
UTIR0540 
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20 IK=IK+K U t IP0550 

C UT1R0560 

30 R0UT(KJS)=-SUM*DINV. UTIR0570 

35 IRLST=IROW UTIR0580 

40 II=Il“IR0W UTIR0590 

RETURN UTIR0600 

50 FORMAT (1H0#10X* *RIN DIAGONAL* »I4» *TS ZFRO») UTIR0610 

END UTIR0620 
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SUBROUTINE WGS (W*IMAXW*lW* JW*DW*U* V> WGS00010 
MODIFIED GRAMM-SCHMIDT ALGORITHM FOR REDUCING WDW(**T) TO UO(j(**T) WGS00020 
FORM WHERE U IS A VECTOR STORED TRIANGULAR MATRIX WITH THE WSS00030 
RESULTING D ELEMENTS STORE ON THE DIAGONAL WGSO0O4O 

WGS00050 


WCIW»Jw) 


IMAXW 

IW 

JW 

DW(JW) 


U ( IW* ( iW+1 ) 72 > 
V(JW) 


INPUT matrix to be reduced to triangular form, 
this matrix is destroyed by the calculation 

IW.LE. IMAXW. and. TW.GT. l 

ROW DIMENSION of W MATRIX 

NO. ROWS OF W MATRIX* DIMENSION OF U 

NO. COLS OF W MATRIX 

VECTOR OF NON-NEGATIVE WEIGHTS FOR THE 
ORTHOGONALIZATION PROCESS. THE D'S ARE UNCHANGED 
RY THE CALCULATION. 

OUTPUT UPPER TRIANGULAR VECTOR STORED OUTPUT 
WORK VECTOR 


(SEE BOOK! 

* FACTORIZATION METHODS FOR DISCRETF SEQUENTIAL ESTIMATION * 
BY G.J.BlERMAN) 

ESTIMATION 

COGNIZANT PERSONS: G.J.BIERMAN/M.W.nFAD (JPL» FEB • 1978 J 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DOUBLE PRECISION SUM*Z*DINV 
DIMENSION W(IMAXWrl)* DW(1), U(l)» V<1) 


WSS00060 
WGS00070 
WGS00080 
WGS00090 
WGS00100 
WGS00110 
WGS00120 
WGS00130 
WGS00140 
WGS00150 
WGS00160 
WGS00170 
WGS00180 
WGS00190 
WGS00200 
WGS00210 
WGS00220 
WG500230 
WGS00240 
WGS00250 
WGS00260 
WGS00270 
WGS00280 
WGS00290 
WGS00300 
WGS00310 
WGS00320 
WGSG0330 
WGSG0340 
WGS00350 
WGS00360 

OU HEPE IS USED AS A WORK VECTQRWGS00370 

WGS00380 

Q Eq, (4 *9) OF BOOK* NEW DW(J) 


Z=Q.D0 
one=i.do 
IWP2=Iw + 2 
DO 100 L=2»IW 
J=IWP2-L 
SUM=Z 

DO 4o K=1*JW 
V(k)=W(J*K) 

U( K^ = DW ( K ) *V (K ) 

40 sum=v(k)*u(k)+sum 

W(J*J)=5UM 
DlNVzSUM 
JM1=J-1 

IF (sUM.GT.Z) GO TO 45 

W(J*.)=0. WHEN DlNV=0 (DINV=NORM(W( . )**2)) 

DO 44 K=1*JM1 

44 W(J*K)=Z 
GO To 100 

45 DO 7q K=1*JM1 

SUM=Z 

DO 50 1=1 *JW 

50 SUM=W(K*I)*U(I>+5UM 

SUM=SUM/DINV 

DIVIDE HERE (IN PLACE OF RECIPROCAL MULTIPLY) TO AVOID 


WGS00390 

WGS00400 

WGS00410 

WGS00420 

WGS00430 

WGS00440 

WGS00450 

WGS00460 

WGS00470 

WGS00480 

WGS00490 

WGS00500 

WGS00510 

WGS00520 
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possible overflow 
DO 60 1=1 »JW 

60 W<K»I>=WCK»I)-SUM*V(I) 

70 W(J»K)=SUM ft Eq,(4.10) 0= BOOK 

100 CONTINUE O STORED IN W(JrK) 

THE LOWER PART OF W IS U TRANSPOSE 

SUM=Z 

DO 105 K=1 * JW 

105 SUM=dW(K)*W(1»K)**2+SUM 

U(1)=SUM 
IJ=1 

DO 110 J=2»IW 
DO 110 I=1»J 
IJrIJ+1 

110 U(lJ)=W(JrI) 

RETURN 

END 


WGS00530 

WGS00540 

WGS00550 

WGS00560 

WGS00570 

WGS00560 

WGS00S90 

WGS00600 

WGS00610 

WGS00620 

WGS00630 

W6S00640 

WGS00650 

WGS00660 

WGS00670 

WGS00680 

WGS00690 

WGS00700 

WGS00710 

WGS00720 

WGS00730 
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